1D Stochastic Inversion of Airborne Time-domain Electromag-netic Data with Realistic Prior and Accounting for the Forward Modeling Error

Airborne electromagnetic surveys may consist of hundreds of thousands of soundings. In most cases, this makes 3D inversions unfeasible even when the subsurface is characterized by a high level of heterogeneity. Instead, approaches based on 1D forwards are routinely used because of their computational efficiency. However, it is relatively easy to fit 3D responses with 1D forward modelling and retrieve apparently well-resolved conductivity models. However, those detailed features may simply be caused by fitting the modelling error connected to the approximate forward. In addition, it is, in practice, difficult to identify this kind of artifacts as the modeling error is correlated. The present study demonstrates how to assess the modelling error introduced by the 1D approximation and how to include this additional piece of information into a probabilistic inversion. Not surprisingly, it turns out that this simple modification provides not only much better reconstructions of the targets but, maybe, more importantly, guarantees a correct estimation of the corresponding reliability.

In the last decades, technological advancements made it possible to move from the simple mineral target detection to more detailed groundwater mapping [21,22] and geological modeling applications [23,24].
Together with the enhanced capabilities of the instrumentation, the strategies to deal with the collected data have also been largely improved. For example, similarly to what happens in seismics (e.g. [25]), the stacking of the recorded transient curves-by means of moving windows whose widths depend on the time-gate [26][27][28]-can be used to preserve high lateral resolution at shallow, while increasing the reliability of the data at depth; this is, indeed, compatible, at the same time: (i) with the higher signal-to-noise at the early-gates, and also (ii) with the larger spatial footprint at the late-gates.

Methodology
The calculation of the response of a physical system is always characterized by some level of modeling error. For example, even in the case of very sophisticated 3D forward modeling tools used to calculate the / responses of an electrical conductivity distribution caused by the excitation induced by a "perfectly" described ATEM system (is it really possible to characterize a priori an acquisition system? [8,13,44]), the parameterization used and the size of the discretized domain might affect the retrieved response. In this respect, every time we use a 1D forward modeling approximation [43] for the inversion of ATEM data (inherently 3D), we introduce some errors that need to be handled; neglecting this source of additional uncertainty would inevitably lead to artifacts, paving the road to successive geological misinterpretations. In the best scenario, the magnitude of the modeling error will be negligible compared to uncertainty in the data, and can be ignored [57]. More often, modeling errors will be significant, and must be taken into account.
Here, we briefly recall a formal framework in which the modeling error is described through a multivariate normal correlated probability distribution, which can be naturally used when the measurement errors are also Gaussian [57 ]. In probabilistic formulations of inverse problems, the goal is to retrieve the posterior probability density function ( | ) measuring the probability of having the model compatible with the measurements . In accordance with Bayes' theorem, ( | ) is proportional to the product between the prior probability density function of the model parameters ( ), and the conditional probability density function ( | ) . Hence, ( | ) ∝ ( | ) ( ) , with ( | ) connecting the measured data and the model parameters, and, in the specific case of a Gaussian noise distribution, can be written as where: (i) is just a normalization factor, (ii) is the forward modeling operator used during the inversion, and (iii) d is related to the data covariance , and, often, by assuming mutually independent data, can be considered equal to ⁄ where the -th component of the vector d is the standard deviation of the -th data component (in the specific case of the ATEM data, [ d ] is the standard deviation of the / value associated with the -th time-gate). If the model parameters are also assumed to follow a Gaussian distribution, then the prior information about the solution can be formalized as follows: ( ) = exp (− 1 2 ( − 0 ) ( − 0 )), in which: (i) is another normalization factor and (ii) the Gaussian is centered on the reference model 0 . In this specific case of Gaussian distributions assumed both for the data noise and the model parameters, we obtain that the maximizer of the probability ( | ) is also the minimizer of the regularized inversion objective functional ‖ ( − ( ))‖ 2 2 + ‖ ( − 0 )‖ 2 2 . In fact, this is a very well-known result (e.g., [36,[58][59]), to some extent, reconciling probabilistic and deterministic approaches; in particular, if − is taken equal to -with being the Tikhonov parameter controlling the relative importance of the regularization term with respect of the data misfit, and a discrete approximation of the spatial derivative-then, the minimization of the objective functional coincides with the standard Occam's inversion [60].
However, the approach discussed in the present research loosens several of these ansätze, and, in the following: i.
we do not restrict ourselves to the Gaussian assumption for the model parameters distribution ( ) as we are going to consider quite general prior distributions defined through the realizations of those distributions and that will be generated via a geologically informed procedure; ii. the = −1 2 ⁄ will not consist uniquely of the component attributable to the noise in the observations, but it will also include a term incorporating the modeling error. In particular, the modeling error will be assumed to be consistent with a Gaussian probability density ( Δ , Δ ) defined by the mean Δ and the covariance Δ . Hence, the ( | ) in Equation 1 will have now the following expression [47] ( | ) = exp (− 1 2 ( + Δ − ( )) Δ Δ ( + Δ − ( ))), with Δ = ( + Δ ) −1 2 ⁄ . By construction, as it will be detailed in what follows, these new terms Δ and Δ will also depend on the prior geological knowledge available about the investigated area.
It is important to stress that despite us demonstrating the effects of taking into account the modeling error within a probabilistic framework, the Δ = ( + Δ ) −1 2 ⁄ can actually be, in a very immediate way, incorporated into a deterministic framework as well. However, the evident advantage of using a stochastic approach is that we can naturally incorporate complex prior information (rather than enforcing simple-e.g., smooth or sharp-constraints) and those pieces of complex information need to be available in any case since they are used for the assessment of the modeling error.
Moreover, given the arbitrariness we can benefit from by defining the prior distribution directly via its samples, in general, we could potentially have the maximum flexibility and even use very powerful strategies as, for example, those based on Multiple-point statistics (MPS) approaches [61]; in those cases, the prior geological information can be formalized by means of the so-called Training Image (TI) that is, basically, representing the conceptual geological model of the expected target subsurface; MPS algorithms can generate samples of the prior that are statistically stationary with respect to the original TI and that can be used as detailed in the following Subsection. In the present research, however, we use other geostatistical strategies to populate the prior (and consistently estimate the associated modeling error). They will be described in the following Section 3 "Results".

Estimation of Gaussian Correlated Modeling Errors
Here, we follow the strategy detailed in [47] to (i) simulate and (ii) quantify modeling errors caused by using a 1D forward as opposed to a full 3D forward for simulating ATEM data.
Firstly, a sample of the underling probability distribution representing the modeling errors is generated. This is done by generating a relatively large sample of Δ realizations of the prior distribution ( ) , of the modeling error associated to the specificth sample of the prior. As will also be discussed in the following, to feed the 1D forward , unidimensional conductivity models have been extracted from the original 3D realizations of the prior in the locations just below the acquisition system positions.
Assuming that the modeling error can be characterized by a multivariate Gaussian distribution ( Δ , Δ ), the mean and covariance can be trivially computed from the samples . Finally, the assessment of the 1D modeling error can be plugged into Equation 2. The number of realizations needed, Δ , and the validity of the Gaussian assumption on the modeling errors will be addressed below.
In the present research, as the best approximation , it has been considered an implementation of the forward modeling discussed in detail in [62]. The simulations have been preceded by a validation phase, in which the 3D forward modeling results have been compared against known solutions assumed to be exact. In our specific case, we performed preliminary tests of the 3D forward against semi-analytical solutions for unidimensional conductivity distributions; the chosen 3D simulation settings have led to mismatches of a few percentage points (generally around 5%). Consequently, as it will be shown later in the paper, the modeling error inherited from the 3D modeling has been assumed negligible with respect to the other noise sources and it has not been further considered in our analysis.

Inversion Strategies
Concerning the deterministic inversion results discussed in the following, as mentioned, we use an Occam's inversion scheme [63] in which each individual / sounding is inverted independently from the adjacent ones and, therefore, the roughness operator in the regularization terms acts only vertically. To be fair, it is true that this specific kind of prior information, formalized by the stabilizer, is not in accordance with the investigated models (as it will be clear in the descriptions of the tests, characterized by abrupt conductivity changes). Nevertheless, the deterministic algorithm retrieves (smooth) models whose 1D responses are in almost perfect agreement with the inverted 3D / data. The level of data fitting used for the inversion of the noise-free data is, for all tests, 0.01%.
The stochastic inversion consists of an independent extended Metropolis algorithm, which is analogous to the standard extended Metropolis algorithm [64], but with an infinitely long step-length; this choice implies that every new model proposal from the prior is a completely independent realization. Accordingly, its implementation consists of a slight modification of the SIPPI toolbox [65,66], making the proposed samples independently drawn from the prior ensemble. The 10 5 realizations of the prior are generated in advance since a subset of them needs to be used for the modeling error assessment; the appropriateness of this choice, and the convergency properties will be discussed in the section "Discussion".

Results
In this section, we perform two synthetic tests of increasing complexity to investigate the effects of including in the inversion process: (i) the proper prior information, and (ii) an estimation of the modeling error. In both cases, we compare our results against a solution provided by a more standard 1D deterministic inversion.

Test 1: 3D Conductivity Distribution with Homogeneous Layers
The first test (Test 1- Figure 1) consists of a conductivity distribution mimicking possible glacial geological settings typical, for example, of Denmark [24,[67][68] and characterized by an intricated network of paleovalleys. Figure 2 compares, for the sounding locations considered in Figure 1, the 3D responses calculated for the entire 3D model against the corresponding 1D measurements that would be obtained by considering exclusively the 1D portion of the original distribution just below the acquisition position. The acquisition specifications are those of a typical VTEM system [11,12] (e.g., each sounding consists of 54 measurements). Not surprisingly, the 1D responses ( Figure 2b) are characterized by abrupt lateral changes associated with the lateral variations of the conductivity model (Figure 2a), whereas, in accordance with the physics, the 3D / data are much smoother ( Figure 2c). Here, we treat the 3D calculated responses as the 'observed data' and invert them by means of both deterministic and probabilistic inversion methods.
In all cases we assume negligible measurement errors during the inversion, as the main purpose is to focus on the modeling error effects. In this respect, it is important to highlight the level that the modeling error can reach: the difference in the 1D and 3D responses (blue line in Figure 2b) can be as big as ~20%, and, very seldom, smaller than several percentage points. Therefore, the size of this mismatch should make accounting for the modeling error unavoidable; on the contrary, as mentioned before, the common practice is to tackle it by discretionally increase the measurement noise. By properly including the modeling error in the inversion, the assumed measurement error could be potentially reduced. Figure 3 shows the result of 1D deterministic inversion-implementing an Occam's regularization strategy. Clearly, the deterministic result fits the observations extremely well (blue line in Figure 3b) and is capable to retrieve the major features of the true model whereas, in some cases, infers deceptive discontinuous reconstructions of the true interfaces (e.g., ~1600 to ~2000 m). It is worth mentioning that the data sensitivity to the model parameters drops below the first reconstructed conductive interface; hence, eventually, all the conductivity variation below that retrieved deep interface would not be considered reliable by skilled interpreters. In this respect, the deterministic result cannot be considered, at least from a practical point of view, much different from the probabilistic result obtained without accounting for the modeling error (and discussed in depth in the following Subsection). Nevertheless, it is true that the retrieved (laterally discontinuous) features (caused by fitting the coherent modeling error-as it will be clear later) might be challenging to be correctly deciphered and might lead to erroneous conclusions.

Stochastic Inversion without Modeling Error Assessment
The stochastic approach allows incorporation of (in principle) arbitrarily complex prior information, as long as realizations from the prior can be generated. In the present research, we generated independent realizations of ( ) representing buried valley structures (similar to those in Figure 1) by means of a Fast Fourier Transform Moving Average (FFT-MA) strategy [69] providing unconditional realizations of a Gaussian random field of the interfaces' locations; in particular, the corresponding mean values are chosen to be uniformly distributed: (i) between 20 and 30 m for the shallowest interface, and (ii) between 65 and 85 m for the deepest, whereas the associated covariance models are characterized by a standard deviations and ranges, respectively, of: (i) 5 m and 100 m, for the first interface, and (ii) 80 m and 500 m for the second.
More details about possible ways of constructing the samples can be found in [65,66] including, for example, some details about the above-mentioned MPS strategies.
Whilst the geometry is varying, the conductivity values of each layer are kept constant realization-by-realization. For clarity, two examples of prior distribution realizations can be seen in Figure 4. It is evident that the parameters defining the realizations of the prior have been selected to be in agreement with our expectations about the geological structures we are dealing with.
It is worth highlighting that, since the inversion is one-dimensional, the actual (1D) models used to feed the inversion algorithm are the individual columns (and the corresponding 1D response) of each 3D realization (e.g., in Figure 4). Therefore, each single 1D prior model can be described by five parameters (the depths to the interfaces and the three conductivities).
Undoubtedly, the stochastic inversion, with such an informative prior, is quite facilitated and is basically reduced to the inference of the locations of the interfaces. As a matter of fact, as Figure 5 shows, by using the Metropolis algorithm [64,70], we can reconstruct quite satisfactorily the first interface at about 25-30 m depth, whereas we generally underestimate the depth of the second layer (except for the last ~1500 m on the right). However, what is most disturbing is that, even in this simple case, the misleading features reconstructed by the stochastic inversion (here, performed without accounting for the modeling error) appear to be almost certain.   Figure 1a taken at the locations of the ATEM soundings (red dots in Figure 1b) -the blue line represents, sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right, in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 1a).

Stochastic Inversion Incorporating the 1D Modeling Error
If instead, we use a subset of the 3D prior samples (Figure 4) to calculate their actual 3D responses and compare them with the data calculated, this time, by means of the 1D forward modeling applied to the 1D conductivity vertical profile in the center of each selected prior realization, we can estimate the appropriate mean Δ and covariance Δ , and incorporate them into the 1D stochastic inversion scheme. The results of the application of this new scheme to Test 1 is shown in Figure 6. Now, the result almost perfectly matches the true model. When the mean depth of the interfaces does not fit the true conductivity change-for example, near the steep lateral variations around 1000 m and 1400 m-the correct location of the boundaries still lies within the uncertainty bands.  (Figure 2a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used (Figure 4), but no modeling error has been taken into account. The orange bands represent the reconstruction uncertainty defined by the mean and standard deviation values (red points and vertical bars) deduced by the retrieved realizations of the posterior ( | ). As for Figure  3a, the black lines show the locations of the interfaces to be reconstructed.  (Figure 1a). In this case, realistic prior is used ( Figure 4) together with the modeling error assessment. As in Figure 5, the orange bands represent the uncertainty (the associated mean and standard deviations are plotted as vertical red bars) calculated by the retrieved realizations of the posterior ( | ). The black lines show the locations of the interfaces to be reconstructed (Figure 2a).
Since we are dealing with 3D conductivity models, it is probably more appropriate to visualize the consequences of the different inversion schemes over several flight lines (Figure 7) across the Test 1 model. If we examine Figures 7b,c, the same conclusions drawn for an individual vertical section are clearly valid also for the other acquisition lines (and, consistently, also for the lateral intra-line resolution-kindly, compare Figures 7e,g).
A direct assessment of the performances of the different inversion schemes in terms of data fitting can be performed by looking at Figure 8b; in general, the relative mismatch between each of the channels of the observed 3D responses and those calculated (via a 3D forward) from the 3D conductivity distribution obtained via the stochastic inversion (with modeling error appraisal) lays within ∓10%. It is worth noticing that the relatively poor data fitting on the right side of Figure 8b is caused by the higher level of heterogeneity of that end of the model, beyond the surveyed area. In fact, Figure 7d shows the rapidly varying morphology for > 2500 m in correspondence of the section considered in Figure 8 (~500 m); that conductivity variation is not in accordance with the inevitable laterally homogeneous extension of the reconstructed solution (Figure 7g). This mismatch   (Figures 7c and 7g). Each of the 54 plotted time-gates (channels) is depicted with a different color.
affects, not surprisingly, mainly the late-gate measurements. The same does not happen on the other end of the section, where the data fitting is particularly good (Figure 8b); in this case, the reason is that, differently from before, the true model continues largely unchanged towards low values, at least for ~500 m (Figure 7d). On the other hand, the 3D responses generated by the conductivity volume retrieved by the 1D deterministic inversion (Figure 8a) demonstrates, once more, that fitting the data with a 1D forward (as in Figure 3b) does not necessarily guarantee that the corresponding 3D calculated data are in agreement with the (3D) observations; indeed, in Figure 8a, the relative misfit between the 3D calculated responses and the measurements ranges approximately between ∓30%.

Test 2: 3D Conductivity Distribution with Heterogeneous Layers
In Test 2, we apply the same strategy to a more elaborated 3D conductivity model. Test 2's model consists of three layers with similar geometries as in Test 1, but, now, characterized by heterogeneous conductivity values (Figure 9). For sake of completeness, analogously to what has been done for Test 1, also for Test 2, we show the 3D responses ( Figure  10c) along the central profile (Figure 10a) of the 3D conductivity model (Figure 9a). In Figure 10b, we display the 1D response as they would result by considering each column of the true conductivity model as independent (for comparison, kindly, see Figure 2b concerning Test 1). Clearly, for Test 2, the importance of taking into account the modeling error is even more evident: the error we introduce when interpreting the 3D model response in terms of 1D data is never below 4% (blue line in Figure 10b). Figure 11a presents the result of Occam's inversion of the 3D data ( Figure 10c). Again, from Figure 11b, it is clear that it is not difficult to perfectly fit the 3D data with 1D responses, even with a model barely capable to get the very major features of the conductivity model to be inferred.
Also for this more complex test, a 1D stochastic inversion can be performed by following the previously discussed Metropolis approach, making use of precalculated samples. The 10 5 realizations of the prior distribution are very similar to those shown in Figure  4 for Test 1; the only difference is that, consistently with the new model to be reconstructed, also for the prior realizations, the conductivity of each layer is allowed to vary. Figure 12 shows the mean map of the posterior distribution obtained without taking into consideration the modeling error; as for Figure 5 (about Test 1), also for Test 2, the retrieved mean model (Figure 12a) does not capture the lateral variation of the top of the deepest layer and smooths the large majority of the incisions out. Perhaps worse than that is the fact that, as shown by the standard deviations of the conductivity of each layer (plotted in Figure 12b), the leveled-out reconstruction is given for (almost) undeniable.   Figure 9a at the locations of the ATEM soundings (red dots in Figure 9b) -the blue line represents sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 9a).   (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used, but no modeling error has been taken into account. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with the mean and standard deviation vertical profile superimposed for each vertical conductivity profile.
On the contrary, if we take into account the modeling error, the mean map derived from the 1D stochastic inversion provides a quite satisfactory reconstruction of the true model with all its complex morphology (Figure 13a), and when the inferred mean conductivity values do not reflect the true distribution, they are associated with high levels of uncertainty (for example at depth ~100 m and ~1300 m in Figure 13b).  (Figure 10c) associated to the conductivity model of Test 2 (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used together with the modeling error assessment. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with, superimposed, the mean and standard deviation vertical profiles.

Discussion
Are 10 5 samples from the prior enough to guarantee the convergence of the stochastic inversions? What about the subset of 500 prior 3D realizations used for the assessment of the modeling error for the geological settings considered? Is the retrieved modeling error Gaussian as it should be to be able to use the proposed inversion scheme? In this section, we try to answer all these legitimate questions.

About the Numerosity of the Prior Samples for the Convergency of the Stochastic Inversion
Since the propositional scheme of the adopted Metropolis approach is based on a finite number of precalculated samples of the prior, it is important to establish if the abundance of those samples is sufficient to guarantee the convergence. In this respect, we run several inversions characterized by an increasing number of prior's samples. From Figures 14f-h, it is evident that, at least for the simple problem of Test 1, convergence is reached with a numerosity of the prior samples of a few thousand. Hence, an abundance of two orders of magnitude higher should reasonably be compatible with our goal. A similar conclusion can be deduced by considering the evolution of the correlation coefficient between the depth of the deepest interface inferred and the true one, as shown in Figure  15.  -h) reports that numerosity. Clearly, here, we are incorporating the information about the modeling error into the inversion of Test 1's data.

Figure 15
Correlation coefficient between the retrieved and the true depth of the deepest layer as a function of the number of precalculated samples of the prior. Clearly, as for Figure 14, we are considering the stochastic inversion with modeling error assessment applied to Test 1's dataset.

About the Numerosity of the Prior's Samples for the Estimation of the Modeling Error
Still considering Test 1's dataset, we try to set up a strategy for the evaluation of the minimum number of realizations of the prior to be considered for an effective assessment of the modeling error. For all the tests performed in the present research, a maximum number of 500 models (similar to those in Figure 4) have been used to calculate the difference and, in turn, the parameters defining the Gaussian probability density ( Δ , Δ ) describing the modeling error. If we consider the behavior of the mean vector Δ calculated for an increasing number of samples number Δ , we can plot the results as in Figure 16a; from that figure, it appears evident that the assessment of Δ reaches significant stability after considering a few hundred realizations of the prior. A similar conclusion can be even more directly deduced by checking Figure 16b, in which the mean of the difference between Δ and its best estimation Regarding the evolution of the covariance matrix Δ , we can, again, study how it changes over the number of prior samples. Figure 17 demonstrates, once more, that for the considered case, Δ =500 should guarantee a reasonable estimation of the modeling error. In particular, Figure 17b shows that the estimation of Δ has already reached convergence after considering ~400 samples.

About the Gaussianity of the Modeling Error
As mentioned before, the proposed approach is based on the assumption that the modeling error is actually Gaussian. In this section, we demonstrate that this working hypothesis is largely met. In this respect, each subpanel of Figure 18 shows one of the 54 histograms (one for each time-gate) of the corresponding component of the vector of the difference between the 3D and 1D responses calculated for the 500 realizations used for the modeling error assessment. In particular, the red lines represent the Gaussian that is better fitting the experimental histograms, whereas the black lines show the Gaussian profile as inferred from the distribution ( Δ , Δ ). It is worth noticing the excellent agreement between the black and red curves, but, more than that, the fact that for the large majority of the time-gates the modeling error is indeed compatible with Gaussian distributions. Surprisingly, the histograms for the late channels (in particular those for going from ~49 to 54-i.e., basically the last row of Figure 18) show some sort of bimodal behavior whose possible justification is not evident. A possible reason might be connected with the finite size of the 3D simulation domain. However, this highly hypothetical guess will need to be investigated and verified. ] of the difference between the 3D and 1D responses for the 500 samples of the prior used for the modeling error estimation. The solid black line is the Gaussian curve better fitting the histogram, whereas the red line is the Gaussian derived by ( Δ , Δ ).

Conclusions
Throughout this paper, two Tests are used to demonstrate the capabilities (and limitations) of stochastic inversion with geologically informed prior distribution to properly reconstruct the complex target conductivity distributions and the associated uncertainties.
In particular, we show how crucial the correct quantification of the modeling error (based on the prior choices) is. In our specific cases, we deal with ATEM data (mimicking the collection of VTEM measurements over geologies recalling glacial sedimentary environments typical, for example, of some regions in Northern Europe) inverted via a 1D stochastic approach with realistic prior and with the assessment of the corresponding modeling error. The conclusion is that, even in the case of stochastic approaches (already improving the results, for example, when compared with 1D deterministic approaches relying merely on the regularization term for including the prior knowledge about the targets), neglecting the fact that we are using a brute 1D approximation, instead of a more sophisticated 3D one, leads to unrealistically low levels of uncertainty in correspondence of those features that might turn out to be just artifacts. On the other hand, the assessment (and utilization) of the modeling error allows a more effective reconstruction of the true models and their associated reliability levels.
This research shows how the initial working hypotheses-concerning the minimum possible number of realizations defining the prior distribution and the numerosity of the smallest subset useful for an effective estimation of the modeling error-can be validated. Moreover, also the ansatz regarding the Gaussianity of the modeling error is largely verified.
Thus, the discussed workflow (presented and tested before on other kinds of data [47]) paves the way to the implementation of 1D probabilistic inversions of ATEM measurements capable to incorporate the complex pieces of geological information available and overcomes many of the difficulties connected with the utilization of efficient 1D approximations.
It is evident that the inclusion of the modeling error does not slow down the already available (and extremely fast) algorithms for 1D inversions (it does not matter if stochastic or deterministic). Hence, this approach will possibly remain useful also when fully 3D stochastic inversions will be practical; in fact, it will not be possible to consider any forward modeling tool perfect and, consequently, accounting for the modeling error will be, most likely, always beneficial.
It is also worth being highlighted that, at least for the investigated Tests, simply a few hundreds of 3D forward simulations are needed to retrieve a robust assessment of Δ and Δ , and, actually, the same estimation for Δ and Δ can be used, in principle, in any survey characterized by similar conditions (i.e. similar prior). Hence, in those cases, the efforts for the calculation of Δ and Δ would impact merely the first survey and clearly would not increase with the size of the survey. On the contrary, a full 3D inversion requires at least a few tens of iterations (i.e. 3D calculations) for each sounding location; this easily results in thousands of expensive 3D forward simulations. From these considerations, it clearly appears how convenient the proposed approach is. On the other hand, it is probably true that severely 3D targets will be, in any case, poorly reconstructed by using 1D approaches; however, the proper inclusion of the modeling error will be always useful to correctly estimate the high uncertainty of the 1D reconstruction of 3D inclusions, whereas not taking into account the modeling error will, most likely, lead to wrong solutions that look (incorrectly) certain.