Linking strike directions to invariant TE and TM impedances of the magnetotelluric impedance tensor

Estimation of the traditional transverse electric (TE) and transverse magnetic (TM) impedances of the magnetotelluric tensor for two-dimensional structures can be decoupled from the estimation of the strike direction with significant implications when dealing with galvanic distortions. Distortion-free data are obtainable by combining a quadratic equation with the phase tensor. In the terminology of Groom–Bailey, the quadratic equation provides amplitudes and phases that are immune to twist, and the phase tensor provides phases immune to both, twist and shear. On the other hand, distortion-free strike directions can be obtained using Bahr's approach or the phase tensor. In principle, this is all that is needed to proceed to a two-dimensional (2D) interpretation. However, the resulting impedances are strike ignorant because they are invariant under coordinate system rotation, and if they are to be related to a geological strike, they must be linked to a particular direction. This is an additional ambiguity to the one of 90° arising in classic strike-determination methods, which must be resolved independently. In this work, we use the distortion model of Groom–Bailey to resolve the ambiguity by bringing back the coupling between impedances and strike in the presence of galvanic distortions. Our approach is a hybrid between existing numerical and analytical methods that reduces the problem to a binary decision, which involves associating the invariant impedances with the correct TE and TM modes. To determine the appropriate association, we present three algorithms. Two of them require optimizing the fit to the data, and the third one requires a comparison of phases. All three keep track of possible crossings of the phase curves providing a clear-cut solution. Synthetic and field data illustrate the performance of the three schemes.


Introduction
The basic unit of the magnetotelluric (MT) method of geophysical prospecting is a 2 × 2 complex tensor per frequency in a given range of frequencies. The measurements are usually taken in a coordinate system that is not the one ultimately used in the final interpretation. In general, the matter of going from one system to the other is far from trivial as explained by Groom and Bahr (1992) in their classical tutorial paper. When it is assumed that the target is a two-dimensional (2D) structure, one of the axes must be parallel to the strike of the structure. It used to be common practice to find the appropriate angle simply by applying the traditional rotation matrix. However, this may lead to erroneous results when the data are polluted by what is known as galvanic distortions (e.g., Bahr 1988). "Rotate at your peril" is the subtitle of a paper by Jones and Groom (1993) where the matter is discussed at length. They use the distortion model of Groom and Bailey (1989) assuming that the local distortions are three-dimensional (3D) but the regional or undistorted impedances are 2D. At present, many MT surveys are interpreted using 3D models (e.g., Ruiz-Aguilar et al. 2020). However, 2D models are still useful in many applications, either to support 3D interpretations (e.g., Martí et al. 2020) or simply because they are the most adequate for the problem at hand (e.g., Comeau et al. 2020).
The Groom-Bailey (GB) approach consists of fitting the data with the response of the model, which includes as unknowns the angle of rotation and two distortion parameters along with the undistorted data, which consist of an anti-diagonal 2 × 2 complex tensor. The algorithm STRIKE based on the extensions made by McNeice and Jones (2001) is the standard for 2D applications. Jones (2012) enumerates four advantages of the approach over other methods. Here, we use a fifth one. It is a rather obvious fact that the undistorted impedances are linked to a strike direction. This coupling is important because if the impedances are to be related to a geological strike they must be linked to a particular direction. The link between the strike estimation obtained from a tensor measurement and the estimate of impedances comes naturally. This fifth advantage of the distortion model

Graphical Abstract
Page 3 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 is a feature missing in other methods for estimating the impedances as discussed in the following paragraphs. For instance, the classical approach of Swift (1967) minimizes the size of the diagonal elements of the impedance tensor to find the optimum strike. Then, the measured tensor is rotated to obtain the new elements in the rotated coordinates. Swift's approach was developed before the community recognized and tried to overcome the effects of the galvanic distortions. Bahr (1988) realized that a strike direction immune to galvanic distortions could be obtained by imposing the same phase on the columns of the impedance tensor. He developed an analytic formula for strike that is immune to galvanic distortions. However, even if one has an undistorted strike and uses it to rotate a distorted impedance tensor, in general, the result cannot be but a distorted tensor because of mode mixing (e.g., Jones and Groom 1993). In contrast, in the GB approach, the fitting process optimizes the strike direction and the undistorted data simultaneously since they are coupled from the beginning. One of the outcomes is that there is no need for further rotations.
Other methods do not need the strike angle to estimate the undistorted impedances but in this case, the impedances are decoupled from the strike. This is possible by the use of invariants under rotation of coordinates. Of the stockpile of invariants of the impedance tensor that have been proposed in the literature, some have a special property when applied to data from a 2D structure: they relate to the traditional TE and TM impedances. One pair of these invariants, max and min , is derived from the phase tensor of Caldwell et al. (2004).
Another pair, represented as complex resistivities ̺ ± , is derived from a quadratic equation (Gómez-Treviño et al. 2014a). We use this last pair in this work. The idea is to restore the strike direction lost when formulating the quadratic equation for the invariants. The issue does not exist if strikes and impedances are treated independently, as in Muñiz et al. (2017). However, when inverting the data, we need to link a strike angle to either ̺ + or ̺ − . This is, of course, impossible on the grounds of how they are computed. The resistivities are invariant under rotation, so they have no specific angle. On the other hand, whether the strike is estimated using Bahr's formula or the phase tensor, it is decoupled from the undistorted impedances. All things considered, the price for going analytical is further uncertain beyond the 90° ambiguity. There is no major difficulty when the modes are easily identifiable, as in a recent application in a marine environment (Montiel-Álvarez et al. 2020). However, if this is not the case, the situation can become helpless. Interpreting the well-known BC87 dataset, much to our dismay, we had to recur to the algorithm STRIKE simply to identify the modes (Gómez-Treviño et al. 2018), although the rest of the analysis was based on the ̺ ± formulae. To identify the modes there are other alternatives. The logical thing to do would be to insert back the ̺ ± analysis into the original GB distortion model. This would come naturally since the resistivities ̺ ± were developed precisely as solutions of this model. In this work, we explore this alternative using several new developments on the subject.

The Groom-Bailey decomposition
Let us begin by assuming that for a given frequency or period, we have a complex impedance tensor Z given as: The impedance tensor is the elementary unit of the MT method, but to make it represent the real Earth some modifications are needed. Small, near surface heterogeneities that are of no interest can severely distort the elements of the tensor. These effects can be modeled using a dimensionless and frequency-independent 2 × 2 real matrix C (e.g., Bahr (1988)). The distorted or measured impedance tensor would be Z m = CZ . The elements of C must be determined or avoided in order to have distortion-free data. The Groom and Bailey (1989) decomposition opts for determining them by first separating those that are determinable from those that are not. The decomposition is represented as The different factors are shown explicitly in Table 1. It is assumed that the undistorted impedances Z 2 are 2D. The rotation matrix R and its transpose R T depend only on one parameter, the strike angle θ . The twist tensor T also depends on one parameter t , as does the shear tensor S with its variable e . Tensor A contains two parameters a and b that are not considered unknowns. They are absorbed as real scaling factors by the undistorted impedances. They account for the well-known static effect on the measured electric field. The full realization that these factors are not determinable from the impedances alone came with the work of Bostick (1984Bostick ( , 1986. The decomposition perfectly acknowledges this fact. Also, something that is not often mentioned is that the factors are applied in the coordinates of the undistorted impedances and not where the measurements were made. However, this is how it should be since the undistorted impedances are the ones that are ultimately interpreted. If data are available for several nearby sites, the static factors can be estimated using existing inverse routines in anti-Occam mode (Gómez-Treviño et al. 2014b, (1) Z = Z xx Z xy Z yx Z yy .
Page 4 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 2018). The problem then reduces to finding the scaled impedances and the distortion parameters twist and shear as well as the strike direction. The solution is obtained by minimizing the fit to the data using the objective function where n T is the number of periods T i , Z mjk is the measured impedance and σ 2 jk is the corresponding variance, and Z cjk is the computed impedance from the distortion model. We use the original objective function defined by Groom and Bailey (1989), although the standard algorithm STRIKE of McNeice and Jones (2001) uses a modified version.
As stated in the Introduction, distortion-free impedances and strikes can be obtained directly from the data without fitting the distortion model. This is done by combining the phase tensor with the quadratic equation. In the following sections, we briefly review both topics concerning how they are used in the present work.

The phase tensor
The phase tensor was introduced by Caldwell et al. (2004) to obtain distortion-free information from distorted impedance data. Separating the impedance tensor Z in terms of its real X and imaginary Y parts such that Z = X + iY , the phase tensor is defined, regardless of dimensionality, as: (3) The distorted impedance Z m can be expressed as Z m = CZ . Separating the product into its real and imaginary parts Z m = CX + iCY , the distorted phase tensor can be written as m = X −1 C −1 CY = . That is, the galvanic distortions that so severely affect the impedance tensor do not affect whatsoever the phase tensor. This property is mainly responsible for its popularity in several applications of the magnetotelluric method. These applications spring from the factorization they make of the phase tensor as: where α − β and R are the strike angle and the rotation matrix, respectively. max and min are the singular values of the phase tensor. Caldwell et al. (2004) developed analytic formulas that provide strikes period by period which in some cases are not consistent showing instability (e.g., Jones 2012). However, it is possible to obtain stable strikes using more than one period at a time as shown by Bravo-Osuna et al. (2021). We use this approach which relies on reframing the phase tensor as: To obtain this equation, solve for the matrix of singular values in Eq. (5) and then assume that θ = α − β is unknown. The angle β is computed using the formula given by Caldwell et al. (2004). The strike angle is obtained by minimizing The summation goes over any desired number of periods T i . Minimizing C L 2 (θ ) overcomes the limiting constraint of assuming that ′ 12 = ′ 21 = 0 for individual periods. The other quantities we use derived from the phase tensor are max and min as given by Bibby et al. (2005) as: The rest of the tensors are defined in terms of one or more parameters that are to be determined in the fitting process Page 5 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 The quadratic equation The GB factorization or decomposition given by Eq.
(2) has some peculiar properties that allow the direct computation of distortion-free impedances, save the static factors a and b . This can be done using the quadratic equation proposed by Gómez-Treviño et al. (2014a).
where ω stands for angular frequency and µ 0 for the magnetic permeability of free space, the distortion-free impedances associated with Z 2 can be computed as The formula is immune to twist and also to the strike angle because ̺ sm and ̺ pm are invariant under rotation. The factor ε −2 within the square root sign in Eq. (12) corrects for the effect of shear. The following section explains how this factor is estimated by combining the phase tensor with the quadratic equation.

Impedance tensor vs. phase tensor
When it comes to directional properties, the GB model and phase tensor are about equally equipped. They hold a similar level of information on the 2-D strike direction. Of the two, we prefer the former because it is less ambiguous. The phase tensor loses part of its uniqueness when the phase curves cross each other. This is because max and min do not keep track of possible crossing of the phase curves. It is possible to recover the original curves as discussed by Booker (2014), but the procedure is far from trivial. By choosing the distortion model, we have to deal with the distortion parameters twist and shear. In the following section, we present three algorithms to link the invariants from the quadratic equation with strikes from the phase tensor. First, we estimate the strike and substitute it into the distortion model. The shear-dependent quadratic solutions are then substituted as the undistorted impedances. A primary scheme considers a straight fit to the data by searching for the optimum twist and shear. Secondly, a simplified version optimizes only for the twist, given that the absolute value of shear can be estimated independently. The third scheme uses the strike and the absolute value of shear estimations, but requires no further optimizations.

Implementation of algorithms
We begin with an example using synthetic data corresponding to site 15 of the distortion-free 2D COPROD2S1 dataset proposed by Varentsov (1998). The data are provided as noiseless apparent resistivities and phases for 12 periods. They are shown in Fig. 1 as continuous and dashed lines, respectively. The colored dots represent the distorted values assuming twist = 20°, shear = 30° and strike = 30°. A small random error of 0.1% of the xy and yx impedances was added to the real and imaginary parts of the distorted data. Notice that there are four distorted curves indicating that the corresponding impedance tensor is 3D, otherwise it would be symmetric. This is obviously an effect of the distortions. The target is to assign directions to recovered undistorted responses that are invariants under rotation of coordinates. In this section, we present three algorithms to achieve this goal. They have in common an independent estimation of strike and differ in how they approach the GB distortion model. Fitting the data with the GB distortion model does not lead to ambiguity of modes, so the basic thing that comes to mind is to mimic the fitting process. To this end, we first estimate the strike using Eq. (7). The penalty function sampled from − 90° to 90° taking into account all 12 periods simultaneously is shown in Fig. 2. We use bootstrapping to estimate the strike and its uncertainty. There are actually 100 penalty functions associated with the same number of realizations, but they are not noticeable in this example. Scattering of the penalty functions will become apparent for more significant errors as seen in later examples. The minima fall as expected at − 60° and 30°. In both cases, the minima are very well defined because of the small random error of 0.1%. We use this strike estimation below for the three algorithms.
Algorithm 1 The first step is to substitute the estimated strikes in Eq. (2). The twist is sampled from − 90° to 90° and the shear from − 44° to 44°. This process represents a joint minimization over the shear and twist. The corresponding values of χ 2 are plotted in Fig. 3. It can be observed that the best fits correspond to Fig. 3b, c. The former indicates that for a strike of 30° ρ + must be associated with ρ yx as is actually the case. The second indicates that for a strike of − 60°, ρ + must be associated with ρ xy , Page 6 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 which is how it should be. Notice that in this case, the minimum falls at a shear of − 30° which is the negative of the assumed shear. It is easy to show that this is how it should be by rotating the right-hand side of Eq. (2) by an angle of − 90°. The operation changes the sign of shear and interchanges the roles of ρ xy and ρ yx . Interestingly, the worst fits shown in Fig. 3a, d have about the correct twist and shear. However, the difference in fitness with respect to the others is so large that they must obviously be discarded.
It is worth noticing that in Fig. 3b the contours of χ 2 are symmetric over the y-axis along the white vertical line. Because of this symmetry the uncertainties in the estimations of twist would tend to be unbiased. On the other hand, because of asymmetry the uncertainties in the shear estimations would need special care. However, to decide the correct association of modes there is no need for detailed analysis of the distorting parameters. The main issue here is the better fit shown in Fig. 3b as compared with that in Fig. 3a. Figure 4 shows explicitly that the choice for a strike of 30° must be ρ + as ρ yx . We show only the apparent resistivities associated with the main components of Z m . The wrong choice presented in Fig. 4a shows an almost antithetical fit of the two modes, while the correct choice of Fig. 4b shows a perfect fit.
Algorithm 2 A second and handier scheme can be implemented considering that the absolute value of shear can be estimated independently, as described in Muñiz et al. (2017). The phases ∅ + and ∅ − from the quadratic equation depend on shear, but max and min from the phase tensor do not. In both instances, the phases are invariant under rotation of coordinates and reduce in 2D to the phases of the TE and TM mode. To compare them, we simply sort ∅ + and ∅ − according to their maximum and minimum values for each period and take the root mean square (RMS) of the difference. Figure 5 shows the RMS of the difference for a range of absolute values of shear from 0° to 44°. It can be observed that the minimum of the RMS clearly recovers the absolute value of shear. Notice that in this case, there is no need to specify the strike as in Algorithm 1. We can now proceed to optimize for twist considering the two possible signs of shear. This is a more straightforward operation that should lead to the same result and in fact, it does. Figure 6a presents the worst fits considering the strike of 30° (Fig. 2) and assuming that ρ + is associated with ρ xy . Neither sign allows for a good fit. On the other hand, as shown in Fig. 6b, assuming for the same strike that ρ + is associated with ρ yx offers a clear-cut picture. The sharp minimum corresponds to the correct positive shear and to the assumption that ρ + is associated with ρ yx , which is the same conclusion reached by Algorithm 1.
Algorithm 3 A third algorithm that still relies on the GB model but does not require fitting the measured data can be implemented as follows. Consider the original GB decomposition and multiply on the left-hand side by R T and on the right by R . The result is: The diagonal elements of the tensor Z ± are zeroes. Z + occupies the place of Z xy and Z − that of Z yx as before. The matrixes T and S are real so their product can be written as: The constants c, d, f and g are real but otherwise, they are arbitrary. Equation (16) can be written as: The elements Z Rij on the left-hand side stand for the rotated elements of the measured or distorted impedances. These can be computed considering the strike has been estimated independently. The right-hand side tells us how twist and shear distort Z + and Z − , which have been corrected for shear through the quadratic equation. The fact that the constants c, d, f and g are real implies that twist and shear affect the magnitudes of Z + and Z − but not their phases. This means that regardless of distortions, the phase of Z + should be the phase of Z Rxy if it corresponds to the chosen strike. If not, it should be associated with the phase of Z Ryx . Figure 7 shows the difference of phases between the two alternatives. There is no doubt that for a strike of 30°, Z + must be associated with Z Ryx . This is because the difference of their phases is nil compared with the choice of associating Z + with Z Rxy . Explicit comparison of the measured and computed values for a strike of 30°. a Assuming that ρ + is ρ xy and b assuming that ρ + is ρ yx . We present only the corresponding components xy and yx. It can be observed that the computed resistivities in b fall in the correct sense. This means that we must identify a strike of 30° with ρ + as ρ yx Page 9 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 At this point, it is worth remarking that the strike does not enter into the computation of the impedances ( Z + and Z − ) because they are invariant under rotation of coordinates. This means that errors in the determination of strike do not transfer into errors in the impedances. This is one of the assets of the proposed algorithms that contrasts with the standard approaches of fitting the data simultaneously for all parameters. Figure 8a presents the penalty functions for the estimation of strike assuming errors of 5% on the original distorted data. They correspond to 100 realizations when using all 12 periods in Eq. (7). The average of the  The other parameter that is needed to apply Eq. (18) is shear. As explained above, this is estimated by comparing the phases ∅ + and ∅ − with max and min . This is done for sampled values of shear from 0º to 44º. The differences for the 100 realizations are shown in Fig. 8b where it can be observed that the estimate of the absolute value of shear is 28.64º, which is close to the true one of 30º, although not as close as the strike estimation. A better estimate can be obtained by increasing the number of realizations. This is important for the accurate computation of Z + and Z − but as illustrated further below, a very accurate value leads to the same binary decision.
The decision is made by comparing the phases of Z + with those of Z Rxy and Z Ryx . Figure 9a presents the differences for the 12 periods. It can be observed that except for two periods, the third and fourth, the differences clearly indicate that Z + must be associated with Z Ryx . This is reflected in the smaller RMS of 2.9º against 29 o for the other option. Although in principle, only the phases for one period are needed to make a decision, it is safer to use many of them. Figure 9b presents the actual Fig. 7 a Differences of the phases of Z Rxy and Z Ryx and ∅ + . b Explicit comparison of phases for a strike of 30° assuming that ∅ + is associated with the phase of Z Ryx  -Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 comparison of ∅ + with the phases of Z Ryx and ∅ − with the phases of Z Rxy . Notice that there is a consistent more significant scattering of the phases of Z Ryx as compared with the estimates of ∅ + . This is because the rotation matrixes in Eq. (16) amplify the random noise while the invariant phases require no knowledge of the strike. Figure 10 explicitly shows the phases of Z + and Z − for the estimated shear of 28.64º as they compare with those corresponding to the actual value of 30º. It can be noticed that the differences are not significant and that most estimates fall each other within their corresponding error bars. In this case, the binary decision, i.e., Fig. 9 a Differences of the phases of Z Rxy and Z Ryx and ∅ + . b Explicit comparison of phases for a strike of 30° assuming that ∅ + is associated with the phase of Z Ryx . Notice the larger scattering of the phases of Z Ryx as compared with the values of ∅ + Fig. 10 Phases ∅ + and ∅ − for the estimated shear of 28.64º as they compare with those corresponding to the true value of 30º Page 12 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 which invariant corresponds to which mode is hardly affected by the accuracy of the shear estimate.

Application to field data
We chose site lit902 of the established BC87 dataset (e.g., Jones 1993) to illustrate the application of the algorithms to field data. The apparent resistivities and phases are shown in Fig. 11. It is possible to observe that they have a 3D character, probably because of galvanic distortions. Clearly the impedance tensor is not symmetric. Figure 12 shows the estimation of the strike using the reframed phase tensor through minimization of the penalty function given by Eq. (7). The estimate is based on using all the available periods and sampling the function from 0° to 90°. Although there is some scatter the minimum is well defined around 62°.
Algorithm 1 Following the steps described above for synthetic data, the strike estimation is substituted into Eq. 2. Again, twist is sampled from − 90° to 90° and shear from − 44° to 44°. The penalty functions χ 2 considering a strike of 61.87° and assuming Z + as Z xy and Z + as Z yx are plotted in Fig. 13. It is observed that the best fit corresponds to Fig. 13a indicating that for a strike of 61.87° Z + must be identified as Z xy . Figure 13b assuming Z + as Z yx still shows a minimum but it is larger than that of Fig. 13a. This result also indicates that for a strike of − 29.13° Z + must be associated to Z yx .
Algorithm 2 This algorithm requires alongside the strike estimation, the knowledge of the absolute value of shear. Recall that the quadratic equation provides phases ∅ + and ∅ − that depend on the square of shear, and that the phases max and min do not depend on shear, so their difference should disappear when the correct value of shear is used. Again, the penalty function, in this case, is the norm of the difference between the phases of the two approaches. The estimated absolute value of shear is 6.9°. The penalty functions are shown in Fig. 14a. Although there is some scattering, the minimum is well defined.
The next step in this algorithm is to determine the twist and the sign of shear. The χ 2 objective functions assuming Z + as Z xy are shown in Fig. 14b for 100 realizations. It is noticeable that unlike the example with  Additional tests which we do not show here indicate that this is due to the noisy data from short and long periods. Removing these noisy data produces single minima. Despite this splitting of the values of χ 2 , it is clear that the negative sign produces the lowest minimum. Figure 14c shows the counterpart assuming Z + as Z yx . The same splitting of minima is present for both positive and negative shear. The ambiguity of whether Z + corresponds to Z xy or Z yx is decided by comparing the minima in Fig. 14b, c. The lowest Fig. 13 The contours represent the misfit to the data: a strike of 61.87º and Z + as Z xy with a smaller misfit. b Strike of 61.87º and Z + as Z yx with a more significant misfit Fig. 14 a Estimation of the absolute value of shear comparing the phases from the quadratic equation and the phase tensor. χ 2 objective functions for the estimation of twist: b assuming Z + as Z xy and c assuming Z + as Z yx . The black lines in b and c represent the functions for positive shear and the red lines for negative value. The lowest χ 2 corresponds to the assumption of Z + as Z xy and a negative shear Page 14 of 16 Arellano-Castro and Gómez-Treviño Earth, Planets and Space (2021) 73:203 χ 2 corresponds to the assumption of Z + as Z xy and a negative shear. This is the same result obtained above using Algorithm 1.
Algorithm 3 This scheme does not require knowing the sign of shear either to optimize for the twist, so we can proceed directly to compare the phases of Z + with those of Z Rxy and Z Ryx as defined by Eq. (18). The corresponding differences are shown in Fig. 15a. For the computation of the RMS of the differences, we have left out some data from short and long periods as shown by the black rectangle. The computed RMS are 2.9º and 17º for (∅ + − ∅ Rxy ) and (∅ + − ∅ Ryx ) , respectively. This indicates that Z + must be identified as Z xy , which is the same result given by the previous two algorithms. In an ideal 2D model one of the curves must necessarily be around zero. In practice because of 3D induction effects this cannot be fully realized for all periods. In this respect our approach has the same limitations as the GB distortion model. Figure 15b explicitly shows the comparison of the corresponding phases for all available periods. In this case, the ambiguity of whether Z + is identified as Z xy or Z yx is resolved by comparing directly the corresponding phases.
We now come to the end of the process, which is the identification of resistivities and phases with a given mode, regardless of whether we compute all distortion parameters or not. Figure 16 resumes the conclusion of the three algorithms. Apparent resistivities ρ + must be identified as ρ xy for a strike of 62°. And the same for the corresponding phases.

Conclusions
In the Groom-Bailey distortion model all variables are interconnected and to fit the data, it is necessary to revert to elaborated numerical methods. Bahr's strike angle, the phases max and min derived from the phase tensor, and the quadratic equation impedances can be thought of as partial solutions for the Groom-Bailey distortion model. The analytical or semi-analytical nature of these solutions is their main asset. However, their full usage can be severely limited because in the process of breaking the problem into sub-problems, the solutions decouple from each other. The quadratic equation and the phases max and min are invariant under rotation, which means they are decoupled from the strike angle. We can obtain the impedances but do not know to which angle they belong, even if we know the angles. This is in addition to the classical ambiguity of 90°. It is important because the association between strike and impedances provides the physical link with the local geology. Bringing back the original distortion model is the only way to reestablish the connection. All three algorithms reduce the problem to a binary decision by inserting the partial calculations into the original decomposition. We presented three algorithms, two still require optimizing the fit to the data, and the third only needs a comparison of phases on behalf of the Groom-Bailey model. In contrast to the phase tensor where max and min lose track of possible