In this study, the method proposed allows to quantitatively evaluate how stably MT responses can be synthesized using forward calculations. This method is simple and easy to implement, although it requires time to run several forward calculations. It has broad utility, as it can be applied to any numerical forward algorithm or 3D structure models. However, even if the forward calculation is confirmed to be very stable using this method, this does not guarantee the accuracy of the forward calculation because the method does not provide information on the bias component of the uncertainty. Nevertheless, it would be useful to determine the partial uncertainty of the forward calculations for practical 3D conductivity structure models.
The two applications of the method showed that the coefficients of variation of the seafloor MT responses synthesized using FS3D (Baba and Seama 2002; Baba et al. 2013) were between 0.1 and 10 percent for the off-diagonal elements, which is not negligible as they are relatively comparable with typical observational errors. There must be a dependence on the numerical algorithm, but there is no guarantee that other algorithms will always give better results than those presented in this study. In other words, one should not believe that the forward calculation is absolutely accurate. More importantly, the applications clearly demonstrated that the coefficients of variation vary significantly depending on the MT impedance element, period, site, structure model, and coordinate system (Fig. 4–9), suggesting that presuming a constant value as a possible uncertainty in forward calculations is not reasonable.
Here, the quantitative evaluation of how well a given 3D structure model explains the observed MT responses by considering the uncertainty of the forward calculations based on the proposed method. The standard deviation provides information on the extent to which the MT response synthesized by one of \(M\) calculations can deviate from the mean. The information on the uncertainty of the mean estimated from \(M\) forward calculations, that is, how much the sample mean can deviate from the population mean, can be given by the standard error defined as:
$${\epsilon }^{\text{s}\text{y}\text{n}}=\frac{\sigma }{\sqrt{M}}. \left(13\right)$$
The standard error decreases with increasing \(M\) and can therefore achieve a more reliable estimate of the mean with a larger \(M\), in exchange for computation time. The value of \(M\) should be selected based on this trade-off relationship. RMS may be calculated for the residual between the observed MT response and the mean of the synthesized MT responses normalized by the value related to the standard errors for both:
$${RMS}_{2}=\sqrt{\frac{1}{2N}\sum _{i=1}^{N}\frac{{\left|{Z}_{i}^{\text{obs}}-{\mu }_{i}\right|}^{2}}{{\left({\epsilon }_{i}^{\text{obs}}\right)}^{2}+{\left({\epsilon }_{i}^{\text{syn}}\right)}^{2}}}. \left(14\right)$$
Using \({RMS}_{2}\), a squared residual is less or more weighted if the mean of the synthesized MT response is less or more reliable as well as if the observed MT response is less or more reliable. The residuals should be evaluated more reasonably using \({RMS}_{2}\) than using \({RMS}_{1}\).
The impact of the normalization proposed above is demonstrated using the observed and synthesized MT responses for the S Atlantic case. The observed MT responses were obtained from Baba et al. (2017a). Figure 10a shows histograms of the residuals normalized by the conventional and proposed methods. For the conventional case, I calculated the residuals between the observed MT responses that are given in the \(({x}_{\text{N}},{y}_{\text{E}})\) coordinate system and the MT responses calculated for the 6th model that the \(x\)-axis is directed in the north direction (Table 1), and normalized them by \({\epsilon }^{\text{obs}}\). The distribution showed a high concentration in the vicinity of zero and some outliers (the minimum and maximum values were \(-189.9\) and \(+135.9\), respectively, far out of the plot range of Fig. 10a) compared to a normal distribution with a mean of zero and a variance of \({RMS}_{1}^{2}\), \(N(0,{RMS}_{1}^{2})\). In the new method, the residuals were calculated for the observed MT responses and the mean of the 10 synthetic MT responses rotated to \(({x}_{\text{N}},{y}_{\text{E}})\) coordinate system, and were normalized by \(\sqrt{{\left({\epsilon }^{\text{obs}}\right)}^{2}+{\left({\epsilon }^{\text{syn}}\right)}^{2}}\). Although the distribution in the vicinity of zero was similar to that in the conventional case, the number of outliers was significantly reduced (the minimum and maximum values were \(-43.1\) and \(+54.1\), respectively). As a result, the distribution was closer in shape to \(N(0,{RMS}_{2}^{2})\).
The uncertainty of the forward calculation varied significantly depending on the MT impedance element, period, and site, as shown in the two applications. The impact of site dependency on the model evaluation is shown in the site-wise RMS (Fig. 10b). \({RMS}_{1}\) for site Tris11 was extremely large compared to that of the other sites. Most outliers in the residuals mentioned above were residuals for Tris11. The large \({RMS}_{1}\) for Tris11 was primarily attributed to the large absolute residual, because the standard error of the observed response was not markedly different from other relatively good data. The MT response at Tris11 varied significantly with period because of the strong local topographic effect, and the synthesized response to the assumed model did not fit the data well. However, the uncertainty of the forward response was also large, especially in the periods around the cusp of the sounding curves (Fig. 7), resulting in \({RMS}_{2}\) for Tris11 being as small as for the other major sites (Fig. 10b). In the use of conventional evaluation, one must consider reducing the RMS for Tris11 to improve the total RMS; however, this is less meaningful considering the uncertainty of the forward calculation.
Although the implementation of the proposed method for inversion is beyond the scope of this study, I proposed some perspectives in inversion analysis. Suppose that a general iterative scheme of regularized inversion minimizes an objective function. A data misfit should be evaluated by considering the standard error of the forward calculations, as discussed above. This means that the data covariance matrix is dependent on the model; thus, special treatment may be necessary in a precise sense. In addition, it is probably not practical in terms of computational time to conduct more forward calculations in each inversion iteration. A compromise may be to conduct forward calculations for the initial model and fix the covariance matrix during inversion, assuming that the standard error of the forward calculations does not change significantly with a small change in the conductivity model. The standard error of the forward calculations can be updated when the change in the model becomes significant and the updated objective function is minimized. I expect that the update would be necessary more frequently in the earlier stage when the change in the model is generally larger. It is important to avoid falling into the local minimum of the objective function because of overfitting to the data beyond the certainty of the forward calculations.
The error floor, which trims the observed error by a threshold value, is frequently applied in practical MT inversion analyses. \({RMS}_{2}\) may be an alternative to \({RMS}_{1}\) with an error floor. One of the basic motivations for introducing an error floor is to avoid overfitting data with unrealistically small error estimates. However, the criterion of “small” is not always evident. In many cases, it seems largely dependent on the experience of the users, without clear physical, statistical, or numerical evidence. Threshold values are rarely provided to each data point separately; rather, a common value or a relation is applied to particular MT impedance elements for all periods and sites, resulting in the loss of the relative importance of each data point, although there may be cases that are objective. Evaluating the residuals with normalization by standard errors for both the observation and forward calculation proposed above will work similarly to applying an error floor in terms of avoiding the overfitting problem. The advantage is that the evidence and meaning are clear for every data point.
The selection of the coordinate system can affect the model obtained by practical inversion analysis, although the selection is arbitrary for a 3D structure. For example, Tietze and Ritter (2013) reported an example and discussed the causes in terms of the data error interrelated between MT impedance tensor elements and sensitivity to the structure. This study suggests that the uncertainty of the forward calculation can also affect the inversion model variations, depending on the coordinate system. The application of the evaluation method proposed in this study showed that the variance of the synthesized MT responses varied with the rotation angle (Figs. 6 and 9). Therefore, overfitting beyond the certainty of the forward calculation can occur differently, depending on the coordinate system. The trajectory analysis for the rotation of the synthesized MT responses would be useful for discussing the appropriate coordinate system in terms of the uncertainty of the forward calculations.
The proposed evaluation method is also useful for investigating the uncertainty depending on the mesh design, although the above discussion focuses on the uncertainty under a given mesh design. One of the supposed cases is to test whether known structures, such as topography, bathymetry, and/or geological setting, are appropriately incorporated into a conductivity model. Here, I present a test for local topographic effect modeling for site Tris11 in the S Atlantic. I modified the local topography model for the two-stage forward calculation (Baba et al., 2013) by applying a finer mesh design, which expresses the relief of the TDC islands more precisely (Fig. 11a), and conducted second-stage modeling for the same 10 coordinate systems as in the previous modeling. Then, the means, standard deviations, and coefficients of variation of the MT responses rotated to the \(({x}_{\text{N}},{y}_{\text{E}})\) coordinate system were calculated. A comparison with those obtained from the forward calculations using the coarser mesh design (Fig. 11b) showed that the standard deviations and coefficients of variation were mostly improved by using the finer mesh design, except for the \({x}_{\text{N}}{x}_{\text{N}}\) element, although the improvement in the periods in which the MT response showed that the cusp was small for all elements. Furthermore, the mean responses for the coarser and finer meshes were found to agree within their standard deviations, indicating that the bias due to the difference in mesh design is less significant.