Reframing the magnetotelluric phase tensor for monitoring applications: improved accuracy and precision in strike determinations

The magnetotelluric method is increasingly being used to monitor electrical resistivity changes in the subsurface. One of the preferred parameters derived from the surface impedance is the strike direction, which is very sensitive to changes in the direction of the subsurface electrical current flow. The preferred method for estimating the strike changes is that provided by the phase tensor because it is immune to galvanic distortions. However, it is also a fact that the associated analytic formula is unstable for noisy data, something that limits its applicability for monitoring purposes, because in general this involves comparison of two or more very similar datasets. One of the issues is that the noise complicates the distribution of estimates between the four quadrants. This can be handled by sending all values to the same quadrant by adding or subtracting the appropriate amount. This is justified by showing that the analytic formula is also a least squares solution. This is equivalent to define penalty functions for the matrix of eigenvalues and then select the minima numerically. Contrary to the analytic formula, this numerical approach can be generalized to compute strikes using windows of any number of periods, thus providing tradeoffs between variance and resolution. The performance of the proposed approach is illustrated by its application to synthetic data and to real data from a monitoring array in the Cerro Prieto geothermal field, México.


Introduction
The evolution of the magnetotelluric (MT) method, from its scalar to its tensorial version, brought about issues regarding the best ways to process the impedance tensor to obtain information about the subsurface. In many instances what is required is to identify the directionality of the electric currents because this is related to lateral changes in electrical resistivity. In the case of two-dimensional (2D) models, this translates into finding the strike angle for which the measured tensor reduces in some manner to a special case. When one of the axes of the coordinate system is parallel to the strike, the diagonal elements of the impedance tensor are zeroes. Swift's (1967) approach to find the appropriate angle is to assume a rotated version of the measured tensor, which in general is not aligned along strike, and then look for the angle that best fits the 2D criterion. The approach leads to an analytical formula. Another analytic formula that converts the elements of the impedance tensor into a strike angle was developed by Bahr (1988) by imposing the condition that the elements of the columns of the impedance tensor must have the same phases in the case of 2D structures. Groom and Bailey (1989) proposed an approach that numerically fits, in a least squares sense, the measured impedances with an appropriate model that includes the strike angle. In this case the solution is not an analytic formula. Finally, there is the phase tensor of Caldwell et al. (2004) from which an analytic formula is derived for the strike angle.
The approaches described above for strike determinations differ from each other in several ways and are also similar in other ways. Let us consider first the similarities between Bahr's and the phase tensor approaches. Both are devised with the explicit purpose in mind of avoiding the effect of galvanic distortions. Also, they impose rigorous criteria in the sense that, strictly speaking, they can be met only by error-free data. A third similarity is that in both cases the approaches lead to analytic formulae. Let us now consider the similarities between the Swift's and Groom-Bailey's approaches. Both impose somewhat relaxed criteria as opposed to the exact requirements in the other two options. Instead of an exact fit a best fit is sought in the least squares sense. In Swift's case, the antidiagonal elements are not forced to be exactly zeroes but only to have their minimum possible amplitude. The Groom-Bailey approach also complies with a least squares criterion.
Presently it is almost a requirement to include phase tensor ellipses as part of the interpretation of any dataset. Usually they are drawn over maps of the study area to reveal directionality of shallow or deep structures depending on the period of interest (e.g. Martí et al. 2020;Comeau et al. 2020). It is not uncommon for the axes of the ellipses to point in inconsistent directions because of random noise. Any improvement in this respect would be certainly welcome. However, where better determinations of strikes are most needed is in the recent application of the magnetotelluric method in monitoring applications. The preferred formula for monitoring purposes is that derived from the phase tensor (e.g. Peacock et al. 2012Peacock et al. , 2013. One of the reasons is that there is no assumption about dimensionality. All other approaches assume 2D structures. Another reason is that it is immune to galvanic distortions, something that Swift's formula is not. However, its drawback is that it is unstable for noisy data (Jones 2012). This limits its applicability for monitoring purposes to very precise data, to ensure a precise estimation of strike. The question then arises of whether we can improve the estimation of strikes beyond the application of the analytic formula. All the analytic formulae provide strike estimates period by period. That is, each estimate and its variance are independent from the others. However, it is possible to link contiguous periods by assuming smooth variations of strike over period and produce a stable profile in the manner of the Occam philosophy (e.g. Constable et al. 1987). Muñiz et al. (2017) explored this path for the phase tensor using as seeds the estimations that are better constrained. Another possibility for stabilization is exemplified by the Groom-Bailey approach as generalized by McNeice and Jones (2001); the estimates can be made period by period, for a given number of periods or for all periods together. The variance of the estimates generally diminishes as the number of periods increases. In a way, this follows the philosophy of Backus and Gilbert (1968) in the sense that the variance can be improved at the cost of resolution. In this paper we adhere to this viewpoint and look for estimates that can be made period by period, for a given number of periods, or for all periods at the same time.

Phase tensor
The history of the magnetotelluric method is plentiful of examples where the problem to be solved consists of avoiding something undesirable. For instance, consider the chaotic variations of the source strength when making telluric measurements. The variations were neutralized by the inclusion of the magnetic field into the telluric method (Cagniard 1953). In turn, the also chaotic polarizations of the source were taken care of by the impedance tensor (Cantwell 1960). The not less chaotic distribution of small, near surface heterogeneities that produce galvanic distortions was neutralized by the phase tensor (Caldwell et al. 2004). The present work also centers on an avoidance. In this case, mitigate the effect of random noise in the estimations of strike using the phase tensor.
Let us start by a brief account of the development of the phase tensor itself. In relation to distortions, small, near surface heterogeneities mean that depth, size, and distance to the electric line must be smaller than the skin depth. All things considered, local electromagnetic induction is small so that only galvanic effects due to electric charges are important. Physically, the source of the charges is the undistorted electric field, so their strength must be proportional to the strength of the source. In short, the measured or distorted electric field E m can be modeled as E m = CE, where C is a distortion matrix and E is the undistorted electric field. In terms of x and y components: The elements of the distortion matrix C are real on behalf of the large skin-depth assumption. The charges are in phase and follow the undistorted electric field. The diagonal elements account for the distortion produced by a component to the same component, while the off diagonal for the orthogonal contributions. To see how the distorted electric fields affect the impedance tensor, consider the undistorted fields expressed in terms of the impedances. (1) where H x and H y are the corresponding components of the magnetic field and Z ij represents the elements of the impedance tensor. The distorted impedance tensor can be written as follows: The real elements of the distortion matrix affect only the amplitude of the impedances, but this is solely for the individual products. The actual linear combination can produce distorted impedances that resemble none of the originals. This explains the importance of the subject and the attempts to neutralize the effects of the distortion matrix. The formula developed by Bahr (1988) for the strike of 2D structures might be considered a precursor to the more general formula derived from the phase tensor. It happens that for the strike angle the product of the distortion matrix and the impedance tensor has a special property: the elements of each column must have the same phase. Imposing this condition on the measured impedance, Bahr (1988) developed his formula for the strike which holds regardless of the distortion matrix. The phase tensor of Caldwell et al. (2004) generalizes this formula and provides other parameters equally immune to distortions. Separating the impedance tensor Z in terms of its real X and imaginary Y parts such that Z = X + iY , the phase tenor is defined, regardless of dimensionality, as Explicitly where det (X) = X 11 X 22 − X 21 X 12 is the determinant of X.
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 have no effect whatsoever on the phase tensor. This property is largely responsible for its popularity in several applications of the magnetotelluric method. One of these applications is determining the strike direction. To this end, Caldwell et al. (2004) developed their analytic formula for the strike θ as: where and The formula for strike provided by the phase tensor reduces to that of Bahr's (1988) when β = 0 . It is in this sense that we mentioned above that the latter is a special case of the former. Our interest here is in the more general formula, particularly in relation to its vulnerability to random noise. Our curiosity arose by analyzing the performance of several formulae for strike as presented by Jones (2012).
To illustrate step by step how the formula for strike is unstable, we prepared a sequence of figures using synthetic data that correspond to site 15 of the COPROD2S1 dataset (Varentsov 1998). Figure 1 shows the elements of the impedance tensor which were distorted from the originals using the decomposition of C given by Groom and Bailey (1989). Table 1 summarizes this decomposition. The 2D undistorted impedance tensor Z 2 is multiplied by the different operators to obtain the measured impedance Z m , which in this case is not actually measured but simulated. We chose the distorting parameters called twist and shear which along with the strike are expressed in degrees. The particular values are twist = 20° and shear = 30°, and an overall rotation of coordinates of − 30°, so that the targeted strike is 30°. In addition, a relatively small random noise with normal distribution was added to the resulting impedances. The noise was computed as 1% of the off-diagonal impedances and added to all the elements. As stated earlier, the object of the exercise is to explore how these small errors amplify when computing the strike angle.
The first level of computations involves using the impedances in the formula for the phase tensor given by Eq. 5. The computed diagonal elements are shown in Fig. 2a, where it can be observed that they are almost error free for the first six periods and are noisier for the last six. The off-diagonal elements shown in Fig. 2b behave very much the same. Notice that for the first six periods the elements are equal, indicating that the tensor is symmetric. In fact, they should be equal for all periods because the data were drawn from a 2D model. The fact that they are not for the last six periods is due entirely to the noise and to the distortions which introduce 3D effects.
We now turn to the ratios that define the angles α and β . The former involves the sum of the off-diagonal elements as numerator and the subtraction of the diagonal elements as denominator. Notice in Fig. 3a that both quantities go through zero at a period of about 10 s. This will have implications when computing the ratio. Figure 3b illustrates the variations of the sum of the diagonal and the subtraction of the off-diagonal elements. This subtraction falls around zero for all periods, but no problem is expected when taking the ratio because for the angle β the subtraction goes into the numerator. Notice that for the last six periods the Fig. 1 The elements of the impedance tensor correspond to site 15 of the COPROD2S1 dataset. The original impedances were distorted and random noise added as indicated. a Real parts. b Imaginary parts. The values were computed as the mean of 30 realizations per period with their corresponding standard deviations. The continuous lines join the error-free values which in most cases are practically identical to the mean values. Black corresponds to X 11 and Y 11 ; blue to X 12 and Y 12 ; green to X 21 and Y 21 ; and red to X 22 and Y 22 subtraction has less scattering than the individual elements shown in Fig. 2b.
The behavior of the ratios themselves is illustrated in Fig. 4a. Both ratios are about uniform for all periods, as they should since they correspond to a uniform strike. The departures from uniformity are due to the effect of noise. The first observation about the ratio for the angle α is the relatively high noise for the period of 10 s and its two neighbors. This is where the diagonal elements 11 and 22 are very similar to each other, likely because of sensing a 1D structure. Both the numerator and denominator are close to zero for these periods and in particular the denominator. Turning now to the longest periods it can be observed that the ratio is very noisy for the last three periods. Again, this is a reflection of the diagonal elements approaching each other, and likely in this case because of a diluting effect of the 2D structure in view of the increasing skin depth. Figure 4a also shows the ratio associated with the angle β . In this case the ratio is a lot less noisy mainly because the denominator is now the sum of the diagonal elements and not its difference as for the angle α . Thus the strike angle θ = α − β owes its sensitivity to noise mainly to the estimation of the α angle. The last step is to compute actual values of the strike and do the statistics. We compute the strikes for a number of realizations and take their arithmetic mean and its standard error for each period. The results are shown in Fig. 4b. We present three scenarios: (1) random noise of 0.01% using 30 realizations, (2) random noise of 1% using 30 realizations, and (3) random noise of 1% using 1000 realizations. The first essentially recovers the true strike because of the small error, the second recovers approximately the true strike only for the periods anticipated by the analysis shown in Fig. 4a, and the third illustrates that increasing the number of realizations improves both accuracy and precision. However, one cannot but notice that for the period of 10 s the strike is underestimated in spite of the small error and the 1000 realizations. This point deserves further attention because it may be the tip of the iceberg of something deeper, in view that the ratios and the trigonometric function tangent are nonlinear combinations of the data and thus not easy to visualize.
It is difficult just by inspection to predict the outcomes of a nonlinear formula. Consider for instance Fig. 4a. There is a bias toward higher values for some of the realizations around the period of 10 s. One would think that increasing the number of realizations would simply fill out the empty spaces. In fact, this is what happens except for a small but significant difference. Figure 5a presents the estimates for the same 1% error and 1000 realizations. It can be observed how some realizations are consistently negative for the period of 10 s. These values seem isolated outliers but actually they are not. Figure 5a also presents the case of 5% error which illustrates that this is a systematic effect associated with how the strike is computed, and that it must be associated with a change in quadrant. At this level of error the results of using 30 or 1000 realizations are comparable as illustrated in Fig. 5b. This unconstrained use of the analytic formula cannot be used for doing statistics because of the bias introduced by the jump in quadrants.
The overall analysis in this section applies only to site 15 of the COPROD2S1 dataset. It does not pretend to be a general proof that one ratio or the other is most sensitive to random noise. We see it as a window into what is behind the unstable character of the formula for strike. In order to improve upon the performance of the formula, we must deal with two issues. One is to neutralize the effect of the noise in relation to the change of quadrants. The other is to be able to use several contiguous periods at the same time which could further decrease the scattering of the estimations.

Reframing the phase tensor
Let us begin by recalling Swift (1967) formula for strike in terms of the elements of the impedance tensor. The first thing is to consider a rotated version of the impedance tensor through a rotation matrix R . That is, The second step is to solve for the strike angle using a least squares criterion. That is, finding θ minimizes the penalty function. This last step acknowledges that the data are not ideal, in the sense that they may have some errors. The penalty function C L 2 (θ) has a minimum at all quadrants, and given the periodicity of the rotation matrix they are 90° apart from each other. The formula for the strike θ = α − β has the same periodicity, but there is no warranty that in the presence of noise the computed values (10) correspond to some sort of minima. Swift (1967) analytic formula is The asterisk stands for the conjugate of a complex number. It is stated in articles and textbooks that this fol- This seems to be the case of Eq. 11 because Sims and Bostick (1969) obtained the same formula without recurring to the least squares criterion. They consider a special combination of the tensor elements and show that the quantity describes an ellipse when rotated. Then they find the angle of the major axis of the ellipse. The result is Eq. 11. In what follows we show that Eq. 6, which is considered an exact, analytic solution for strike, can also be considered as a least squares solution.
The expressions for α and β in Eq. 6 were obtained originally by Bibby (1977) decomposing the direct-current resistivity tensor into symmetric and skew-symmetric parts. The angles are used by Caldwell et al. (2004) to factorize 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. Notice that the rotation on the left is for the difference of the angles and that that on the right is for the sum of them. To make this product symmetric multiply on the right by R(2β) so that, We now solve for the matrix of eigenvalues as: Assuming that the strike θ = α − β is an unknown, the matrix of singular or principal values is no longer (14)

Fig. 5
The analysis of the formula for strike indicates that it is biased with respect to random noise: a the individual strikes for each realization are plotted for two levels of error. The green lines correspond to 1% error and 1000 realizations. The yellow lines correspond to 5% error and 1000 realizations. The data were distorted assuming twist = 20° and shear = 30°. The dashed line represents the true strike as computed with 0.01% error. b The black line represents the estimates of strikes using 5% error and 30 realizations. The red line represents the estimates using 5% error and 1000 realizations, and the green dashed line corresponds to 0.01% and 30 realizations diagonal unless it corresponds to the appropriate strike. In general, Posed in these terms, the analytic formula for strike obtained by Caldwell et al. (2004) would correspond to setting ′ 12 = ′ 21 = 0 . Using Swift's criterion this would correspond to minimize.
In other words, for every period T i , we would look for the angle θ that minimizes the size of the anti-diagonal elements of the matrix defined by Eq. (14). One way to do it is to proceed as for Eq. 10 and obtain a formula for the least squares solution. This is done in Appendix with the result that the formula is the same as the well-known formula for strike given by Eq. 6. The other way to proceed is to compute the angle β using Eq. 8, and sample the penalty function for different values of θ and select numerically the angle where the function has its minima. This we call reframing the phase tensor for strike determinations.
Evaluating the analytic formula or locating the minimum of the penalty function should provide the same strike. This is in fact so but the way they do it is entirely different. This is illustrated in Fig. 6 for two of the realizations shown in Fig. 5, but using impedances from a single period. In other words, the data are a single impedance tensor altered by random error in two different ways. The target is a strike of 30°. We plot the penalty function from − 90° to 90° and the strikes provided by the analytic function. It can be observed that the two values provided by the formula fall exactly at one of the minima of the penalty function in each case. The important thing here is that even if the estimate falls outside the target quadrant, it falls exactly at another minimum. There is a minimum at each quadrant distant 90° from the contiguous one.
To do proper statistics using either the analytic formula or the minima of the penalty function, we need to place all estimates in the same quadrant. In the case of the formula this can be done by identifying the quadrant of a given estimate and adding or subtracting the proper amount to send it to the desired quadrant. This is illustrated in Fig. 7a for seven realizations where all are sent to the first quadrant. This we call a constrained application of the analytic formula, as opposed to unconstrained when simply accepting the computed values. In the case of the penalty function the moving around of the estimates is not necessary. It is sufficient to limit the sampling to any of the quadrants, in this case the first. Figure 7b compares the performance of both approaches using 30 realizations. It can be observed that the difference between the corresponding estimates and their Fig. 6 Typical migration of quadrants due to errors in the data when using the analytic formula. The continuous lines represent the penalty functions C L2 (θ ) for a single period for two realizations. The circles correspond to the estimates using the analytic formula. The black circle falls on the first quadrant exactly at the minimum of the penalty function and the same happens for the red circle at the fourth quadrant. The dashed lines indicate that the periodicity of the penalty functions is 90° variances are not significant. In fact, they should be identical but are not because they were not computed using exactly the same distributions of 30 samples.
The numerical approach does not add much to the constrained application of the formula, except that it saves the rearrangement of the estimates into a given quadrant. However, it opens the door to the possibility of using in the estimation of strike more than one period, something that the analytic formula cannot handle by itself. To use more than one period we minimize . 7 a Penalty functions and concentration into quadrant I of the strike estimated using the analytic formula. The continuous lines represent the penalty functions for one period and seven realizations. The red circles represent the strikes from the analytic formula. b The black curve represents the estimates using least squares and the red curve using the analytic formula. The green dashed line corresponds to the true strike Bravo-Osuna et al. Earth, Planets and Space (2021) 73:34 The summation is over any number of periods, depending on how wide we desire a windows to be. Figure 8 illustrates how the performance of the phase tensor improves when using windows of several periods. In general, if n T is the total number of periods and n p is the number of periods of the window, there will be n T − n p + 1 windows over which the strike can be estimated. Figure 8a illustrates the case for n T = 12 and np = 2 so that there are 11 windows. The estimates are plotted at the geometric mean of the first and last period. When n p = 1 there are 12 windows, just as when applying the analytic formula. This we call windows of one period. It can be observed in Fig. 8b that wider windows perform even better, Fig. 8 a Least squares estimates assuming windows of two periods. The green dashed line represents the true strike. b Least squares estimates assuming windows of two periods. The green dashed line represents the true strike both in accuracy and precision. We present an example of uniform windows but the application is open for combination of sizes, either fixed beforehand or using nested estimations controlled by their variance. Figure 8 also illustrates that the estimates using windows of several periods are not arithmetic averages of the values provided by the analytic formula shown in Fig. 7b. The joint optimization of several periods performs better than individual estimations or combinations of them.
Another improvement that can be easily implemented is the change of norm in Eq. 17. Analytically, switching from a L 2 to a L 1 norm is a much elaborated mathematical process. However, numerically this is a trivial change. Instead of minimizing the penalty function given in Eq. 17 we can minimize In the presence of outliers it is well known the property of the L 1 norm of providing better estimates, because the L 2 norm by its very nature amplifies the effect of outliers.
At this point it is worth recalling that even though we are applying the same criterion to impedances and to elements of the phase tensor, there are qualitative differences between the two types of data. Swift's formula operates directly on the elements of the impedance tensor. On the other hand, the condition imposed on the phase tensor is applied to combinations of ratios of the impedances. Figure 9a illustrates that working with impedances produces strikes that are more stable than those obtained with the elements of the phase tensor. We used in this case data with no galvanic distortions to isolate the precision aspect of the comparison. Figure 9b illustrates the performance when the data are affected by galvanic distortions. It can be observed that while the size of the errors is still somewhat smaller for impedances, their mean values are way off the target. On the other hand, the strikes from the phase tensor are more accurate but have less precision. This loss of precision seems to be the penalty for being immune to galvanic distortions.
We now turn to explore the recovery of a profile that varies with period. We use this same profile later when comparing recoveries that are very close to each other. The target strikes are 20, 30, and 40°, assuming errorfree data and with small to moderate errors. Figure 10a illustrates the behavior of the penalty function for the error-free case. The minima occur exactly at the 20, 30, and 40° targets, as they should. Although the data have no errors, inspecting the penalty function it is possible to predict which strike value will be the least well determined when including errors. This will be the strike of 20° because it has the flattest penalty function. Then follows 40° and the best determined will be 30° because it has the sharpest penalty functions. This is confirmed by Fig. 10b which corresponds to data assuming a relatively small error of 1% in the impedances. We present penalty functions for only 5 realizations per period to avoid crowding, but somehow these are sufficient to corroborate the predictions mentioned before. It can be observed that the wider scattering of the minima occurs around 20°. Then it is followed by those for 40° and finally by the less scattered around 30°. This is the same sequence anticipated by the relative flatness of the error-free penalty functions. This is in turn reflected in the size of the error bars shown in Fig. 11. The largest errors are for the strike of 20° followed by 40°, and 30° being the best constrained strike. In this case we used 30 realizations.

Application to differential synthetic data
We now turn to the target of detecting changes using the profile discussed above. We place as target a 1° change in strike, using the same site 15 of the COPROD2S1 synthetic dataset, distorted with a twist of 20° and a shear of 30°. Random noise of 5% is added to the distorted impedances. We begin with the unconstrained analytic formula. It can be observed in Fig. 12a that only for the strike of 30° we have the possibility of hitting the target of 1° difference. The other two have too large variances and most estimates are way off the desired targets. Of course, we know that this is not the proper way to use the formula. Figure 12b illustrates the performance of the reframed version using windows of one period. It can be observed that for the group of 20° the estimations are almost as bad as those provided by the unconstrained analytic formula. For the group of 30° the estimations are comparable for the two approaches. This means that for these four periods there are no shifts of quadrants in the analytical solutions. The shifts return for the last four periods. We include this comparison between the unconstrained analytic formula and estimates using windows of one period just for the sake of completeness. In fact, as discussed earlier, the analytic formula when constrained to a single quadrant produces the same estimate as the windows of one period. It is when using wider windows that the variance can be actually reduced.
We present results using windows of even number of periods from four to ten. It is difficult to predict for a nonuniform profile the effect of jumps of strikes when using windows of several periods because of possible incompatibilities. We do not know how the least squares criterion can handle very different values at the same time, and what the effect would be on the variances, which are crucial for detecting changes. The estimations using windows of four and six periods are shown in Fig. 13a. We will only point out that the estimations improve for most of the windows except for the first and possibly for the second. Further improvements are shown Fig. 9 Relative performance of the least squares approach as applied to the impedance tensor and the phase tensor. a With no galvanic distortions. The green dashed line represents the true strike. The black continuous line represents the estimates using the impedance tensor. The red line corresponds to the estimates using the analytic formula derived from the phase tensor. b With galvanic distortions. The green dashed line represents the true strike. The black continuous line represents the estimates using the impedance tensor. The red line corresponds to the estimates using the analytic formula derived from the phase tensor Bravo-Osuna et al. Earth, Planets and Space (2021) 73:34 in Fig. 13b for windows of eight and ten periods. Considering these sequence of results, the overall lesson of the exercise is that we should experiment with windows of different number of periods to obtain changes that are statistically significant.

Application to field data
To illustrate the performance of the approach when using field data, we selected a site from a monitoring network installed around the Cerro Prieto geothermal field, located in northern Baja California, México. Beginning Details about the Pb-PbCl 2 electrodes can be found in Petiau (2000) and Booker and Burd (2006) and about the BF-4-type coils in the operation manual of the MT-1 system (Electromagnetic Instruments Inc. 1996). The sampling frequency was 18.18 Hz. Only two stations were equipped with induction coils but all were synchronized for later processing of the time series using remote reference techniques (Gamble et al. 1979). Specifically, we used the RRRMT8 algorithm developed by Chave et al. (1987) and Chave and Thomson (1989). More details about the network are described by Cortés-Arroyo (2018) and Cortés-Arroyo et al. (2018). We chose to compare 2-week data from 2017 with also 2-week data from 2018 for a station labeled E-4. The comparison is shown in Fig. 14 using windows of one period. In all cases we assumed a 5% error in the data to compute the errors in the strikes. It can be observed that the estimations are very stable for periods between 10 and 100 s. In fact it can be stated that the strike did not change from 1 year to the other in that range of periods. However, for shorter and longer periods the story is very different because of the large oscillations in the strike estimates for both years.
The strikes are a lot smoother when we use windows that include several periods. There is a total of 36 periods, so choosing windows of 6, 12, and 18 would provide 31, 25, and 19 estimates, respectively. The resulting strikes for these windows are shown in Fig. 15. The estimates for the shortest window already show a definite trend for the difference between the two years as illustrated in Fig. 15a. This happens for periods shorter than 10 s, which for the area corresponds to depths of penetration of less than 3-4 km. This would locate the changes within the depth range of the Cerro Prieto geothermal field, with no appreciable change for the deeper and more regional background. It is intriguing that the curves cross each other, implying that the strike went a few degrees to one side for the deeps and to the opposite side for the shallows. This trend is reproduced for the wider windows of 12 and 18 periods as shown in Fig. 15b, c, respectively. However, notice that the negative changes tend to disappear for the widest window, as indicated by the lowering of the red curve of 2017 toward the blue of 2018. Notice also the stability and consistency of the strikes for the longer periods, implying a stable background with a consistent strike for the 2 years. Fig. 11 The estimates of strike correspond to the penalty functions shown in Fig. 10b, except that in this case we used 30 realizations. See the main text for a discussion on how the dispersion of the penalty function relates to the accuracy and precision of the final estimates Fig. 12 Recovery of strikes for two profiles that differ by 1° over all periods: a unconstrained application of the analytic formula. The black line corresponds to the estimates for the base profile of 20, 30, and 40°. The red line corresponds to the estimates of the perturbed profile of 21, 31, and 41°. b Using the reframed version for windows of one period. The black line corresponds to the estimates for the base profile of 20, 30, and 40°. The red line corresponds to the estimates of the perturbed profile of 21, 31, and 41°. In both cases the dashed lines represents the corresponding true values Bravo-Osuna et al. Earth, Planets and Space (2021) 73:34

Conclusion
The unstable character of the formula for strike derived from the phase tensor is not intrinsic to the phase tensor. It is possible to keep the assets of the phase tensor in regard to dimensionality and to its immunity to galvanic distortions by reframing how the strike is estimated. One of the problems to overcome is the effect of random noise which may change the quadrants of the estimated strike. This can be addressed by sending the estimates from any quadrant to a chosen one, which we call the constrained analytical formula. This is justified because the analytic formula for strikes is also a least squares solution. Regardless of the quadrant any estimate falls on a minimum of the least squares penalty function. This means that choosing the minima of the penalty function over any quadrant is equivalent to sending to the same quadrant the estimates of the analytic formula. Using the minima of the penalty function allows the estimation to be made using several periods at the time. This estimation over windows of several periods, something intrinsically impossible using the analytic formula, allows the recovery of strikes with lower variances. Thus, the least squares approach certainly brings something new to the estimation of strike. Regarding uncertainties we did not explore beyond the standard deviation assuming normal distributions.
The subject is open for future work using more elaborated methods as stochastic, Bayesian, or Markov chains approach for sampling and estimating parameters and their uncertainties. Another extension would be the implementation of nested algorithms for either imposing continuity of strikes over period, or getting the optimal combination of window sizes for a given variance. (22) tan 4θ = 2 tan 2α 1 − tan 2 2α .