Probabilistic Petrophysical Reconstruction of Danta's Alpine Peatland via Electromagnetic Induction Data

Peatlands are fundamental deposits of organic carbon. Thus, their protection is of crucial importance to avoid emissions from their degradation. Peat is a mixture of organic soil that originates from the accumulation of wetland plants under continuous or cyclical anaerobic conditions for long periods. Hence, a precise quantification of peat deposits is extremely important; for that, remote‐ and proximal‐sensing techniques are excellent candidates. Unfortunately, remote‐sensing can provide information only on the few shallowest centimeters, whereas peatlands often extend to several meters in depth. In addition, peatlands are usually characterized by difficult (flooded) terrains. So, frequency‐domain electromagnetic instruments, as they are compact and contactless, seem to be the ideal solution for the quantitative assessment of the extension and geometry of peatlands. Generally, electromagnetic methods are used to infer the electrical resistivity of the subsurface. In turn, the resistivity distribution can, in principle, be interpreted to infer the morphology of the peatland. Here, to some extent, we show how to shortcut the process and include the expectation and uncertainty regarding the peat resistivity directly into a probabilistic inversion workflow. The present approach allows for retrieving what really matters: the spatial distribution of the probability of peat occurrence, rather than the mere electrical resistivity. To evaluate the efficiency and effectiveness of the proposed probabilistic approach, we compare the outcomes against the more traditional deterministic fully nonlinear (Occam's) inversion and against some boreholes available in the investigated area.


Introduction
Peatlands are important carbon reservoirs, and some of them are located in mountain settings for which efficient and reliable methodologies for the detection of their geometry, bottom morphology, and volumes are still lacking (Pezdir et al., 2021;Silvestri, Christensen, et al., 2019;Trappe & Kneisel, 2019).Here, we explore the use of ground-based frequency-domain electromagnetic induction (FDEM) technology to quantify peat thickness and extension in a bog located in the Italian Dolomites in which peat overlays electrically conductive clay substrate.Electrical resistivity is the main physical property retrieved via the inversion of the FDEM measurements (Guillemoteau et al., 2016(Guillemoteau et al., , 2019)).Unfortunately, per se, this physical property of the investigated medium is of little importance to map peatlands.Instead, a direct estimation of the probability of the occurrence of peat or clay would be preferable.In this paper, such probability is inferred by means of a probabilistic petrophysical inversion of the geophysical measurements (Bosch et al., 2010;Grana, Russell, & Mukerji, 2022;Guo et al., 2023).
Probably, the first attempts to include petrophysical relationships into probabilistic inversion frameworks date back to the beginning of the twenty-first century, when such schemes were proposed for the uncertainty assessment of the spatial distribution of physical properties in seismic reservoir characterizations (Connolly & Kemper, 2007;Malinverno & Briggs, 2004).Since then, they gained increasing interest over the years (Hansen et al., 2012;Minsley et al., 2021;Zunino et al., 2015).In fact-in addition to directly retrieving the crucial lithological information, rather than merely the physical property-probabilistic petrophysical inversion has the intrinsic advantage of providing an estimation of the uncertainty associated with the results, which is often missing in traditional predictions (e.g., of the electrical resistivity) via deterministic inversions (McLachlan et al., 2021a).Moreover, the proposed probabilistic approach allows for the formalization of quite complex prior information (Bai et al., 2021).This is a further advantage with respect to more standard deterministic schemes in which the ill-posedness of the problem is tackled via the introduction of a regularizing term, which is often capable of formalizing/enforcing (too) simple information/constraints related to the presence of smooth/sharp transitions between different lithologies (Klose et al., 2022(Klose et al., , 2023)).
In the last decades, several ground geophysical methods have been effective in detecting and characterizing peat deposits, both laterally and in-depth.Ground penetrating radar (GPR) has been proven reliable and precise in inferring the deposit thickness (Comas et al., 2015;Warner et al., 1990).Electrical resistivity tomography has also been found successful (Elijah et al., 2012), especially, when used to reconstruct the variations of chargeability associated with the presence of peat (Comas & Slater, 2004;Slater & Reeve, 2002).However, in general, the above-mentioned geophysical methodologies can be applied to relatively small areas and are extremely difficult to perform in flooded environments like wetlands/peatlands (Zaru et al., 2023).In this respect, FDEM presents several advantages as it does not require any direct coupling with the surface.In fact, several FDEM systems are specifically designed to be operated while mounted on helicopters (Karshakov et al., 2017;Yin & Hodges, 2007) and, probably, soon, on drones (Mitsuhata et al., 2022).
Due to measurement and modeling errors, and because of the band-limited observations and finite amount of data points (Qiu et al., 2020), retrieving the spatial distribution of electrical resistivity from the FDEM observations requires solving an ill-posed problem (Tarantola, 2005).In both deterministic and probabilistic inversion approaches, the ill-posedness is addressed by making use of the available prior information.In the deterministic approaches, such prior information is formalized via the regularization term (Zhdanov, 2002).Traditionally, the L2-norm of the spatial gradient of the model parameters is the most common choice for the regularization term (implementing the so-called Occam's inversion- Constable et al., 1987).Hence, in the case of the L2-norm of the gradient, the unique and stable solution of the regularized inversion is selected to be the smoothest possible model that also fits the data within the noise level assumed for the measurements.In fact, the inverse problem is reduced to the minimization of an objective functional consisting of the weighted sum of the regularization term (the regularizer or stabilizer) and the data misfit.The weight between the two terms (generally, indicated with the letter λ and sometimes named Tikhonov's parameter) is selected to ensure adequate matching between observed and calculated responses while minimizing the regularizer (Zhdanov et al., 2006).In many cases, the geological targets to reconstruct are expected not to be smooth; in such cases, the L2-norm of the model gradient should not be used, and more appropriate information (i.e., more consistent with the prior geological expectations) must be included in the objective functional.In this respect, for example, recently, several regularizers based on the minimum support of the model are gaining popularity (Guillemoteau et al., 2022;Rossi et al., 2022;Vignoli et al., 2015).In the case of the minimum gradient support regularization, the inversion algorithm searches for the model with the minimum number of significant variations of the physical parameters (and not for the model with the minimum variation).Thus, the minimum gradient support favors the retrieval of sharp boundaries between the geological structures (Vignoli et al., 2021).
However, the regularization terms need to meet some (mathematical) requirements: not all functional can be regularizers able to transform an ill-posed problem into a conditionally well-posed one (Portniaguine & Zhdanov, 1999).Furthermore, formalizing complex (and realistic) prior information through a set of functionals might be troublesome (if not impossible).In this respect, probabilistic approaches have indubitable advantages: prior information can be described and incorporated into the inversion via the prior distribution, which, in turn, can be defined by its samples (Hansen, 2021a;Hansen & Minsley, 2019).This leads to the possibility of creating the samples in accordance with quite arbitrary expectations about the investigated geology (Bai et al., 2021;Zaru et al., 2023), without the limitations connected to the necessity of formalizing that information in terms of a functional to be minimized.
Another crucial advantage of probabilistic approaches consists in their capability to provide an estimation of the result uncertainty in a very natural way.The posterior probability distribution is the solution of probabilistic inverse problems, and it is typically represented as an ensemble of models (realizations of the posterior probability distribution) conditioned on the data.This is different from the solution of deterministic inversions which consists of the unique model minimizing the objective functional.Hence, in the probabilistic framework, for example, from the final ensemble of models, the uncertainty of a specific feature can be readily extracted by simply counting the relative occurrence of that feature (Hansen, 2021b;Høyer et al., 2017;Minsley, 2011;Mosegaard & Sambridge, 2002;Mosegaard & Tarantola, 1995;Tarantola, 2005).
The above-mentioned inversion strategies are generally used to image the geophysical property distribution (in most of the FDEM methods, that is limited to the electrical resistivity of the subsurface).However, the estimation of petrophysical properties-such as porosity, lithology, and fluid saturations-can be formulated as an inverse problem.In this case, they are usually referred to as petrophysical inversions (in seismic exploration, they are sometimes also known as rock-physics inversion-e.g., Grana, Parsekian, et al., 2022).In the deterministic frameworks, a physical relationship mapping the rock properties to the geophysical parameters is often necessary (Cao et al., 2023;Foged et al., 2014;Mastrocicco et al., 2010).Instead, in the probabilistic context, like the one used in the present research, the connection between lithology and physical properties is formulated in statistical terms and inherently incorporates the associated uncertainties (here, the uncertainty is not the one characterizing the geophysical measurements, but, rather the one associated with the petrophysical link) (Bai, 2022;Grana, 2016Grana, , 2020;;Grana & Della Rossa, 2010;Gulbrandsen et al., 2017;Madsen et al., 2023).
Thus, in the present work, we test a probabilistic petrophysical inversion on an FDEM data set originally collected for the characterization of an Alpine peatland.We assess the performance of the proposed inversion scheme at the locations of the several available boreholes (not used to constrain the inversion) and by comparing the obtained results against more traditional deterministic Occam's inversions.

Materials and Methods
The forward modeling tools used in electromagnetic inverse problems are, nowadays, fully nonlinear and, despite 3D forward options could be available (Cox et al., 2012;Oldenburg et al., 2020), the use of 1D approximations is still the standard in the industry and in the academia (Bai et al., 2020;Dzikunoo et al., 2020;Hansen & Finlay, 2022;Ley-Cooper et al., 2015;Viezzoli et al., 2008).The reason is mainly related to the computational costs of 3D modeling (Cai et al., 2022;Guillemoteau et al., 2017;Kara & Farquharson, 2023).In the present research, we make use of the fully non-linear 1D forward modeling discussed, for example, in Elwaseif et al. (2017a).Given a sequence of layers with known electrical, homogeneous properties-represented by the model m-the 1D forward modeling F(m) maps the resistivity (the sole electrical property we are investigating) distribution into the observed FDEM response collected at the surface (described by d): (1) The goal of probabilistic formulations is to infer the posterior probability density (PDF) function, p(m|d obs ), measuring the probability of the model m being compatible with the observations d obs .According to Bayes' theorem, p(m|d obs ) is proportional to the product of the prior PDF of the model parameters p(m) and the conditional PDF p(d obs |m).Thus, p(m|d obs ) ∝ p(d obs |m)p(m), with p(d obs |m) connecting the measured data and the model parameters.In the specific case of normally distributed noise in the measurements, the conditional probability can be written as Earth and Space Science 10.1029/2023EA003457 ZARU ET AL.
where: (a) k d is simply a normalization factor and (b) C d is a matrix taking into account the correlated Gaussian data noise.Here, we assume the noise to be uncertain and, consistently, consider diag(σ d ) 2 = C 1 d in which the ith component of the vector σ d is the standard deviation of the ith data component.In our specific case, [σ d ] i is the standard deviation of the quadrature component of the ratio between the induced and the primary fields associated with the ith frequency.In fact, during FDEM surveys, the measurements consist of the real and imaginary parts of the ratio between the secondary (interacting with the subsurface) and the primary field (generated by the source).Here, we limit our analysis to the imaginary part (the quadrature component) since it is the one mostly sensitive to the electrical resistivity, whereas the real (in-phase) component is mainly related to the magnetic susceptibility of the investigated medium (Hendrickx et al., 2002).
Assuming a Gaussian distribution also for the model parameters, the prior information about the model m can be formalized as follows: with: (a) k m being another normalization factor and (b) the Gaussian being centered on the reference model m 0 .
When Gaussian distributions are assumed both for the data noise (with Gaussian uncertainties described by the covariance operator C d ) and the model parameters (with prior information on m described by the a priori reference model m 0 with the associated uncertainty represented by the covariance operator C m ), the maximization of the probability p(m|d obs ) corresponds to the minimization of the regularized inversion objective functional . This is a well-known result (e.g., Tarantola, 1987;Tarantola & Valette, 1982;Zhdanov, 2002) that bridges the gap between probabilistic and deterministic approaches.For example, if the inverse of the model covariance C 1 m is taken equal to λ 2 L T L-with the positive scalar λ weighting the relative importance of the regularizer and the data misfit, and L being a discrete approximation of the spatial derivative-then, the minimization of the objective functional implements the standard Occam's inversion (Constable et al., 1987).In this specific case, the objective functional P (λ) (m, d obs ) to be minimized during the deterministic inversion is: On one hand, the strategy we are proposing here relaxes one of these working hypotheses; in particular, in the following, we do not make any assumptions on the Gaussianity of the model parameters distribution p(m) as we consider quite general prior distributions defined via their realizations, which are generated through a geologically informed procedure.
On the other hand, we keep considering the noise contaminating the observed quadrature components to be frequency-by-frequency, and location-by-location, independent, even though it is clearly not (Liu et al., 2023;Minsley et al., 2012;Mitsuhata et al., 2001).Consistently with this assumption, we still approximate C 1/ 2 d with a diagonal matrix, whereas, for example, with the availability of a 3D FDEM forward modeling, we could have included the noise component attributable to the modeling error associated with the 1D modeling approximation.Bai et al., 2021 showed, for the time-domain electromagnetic case (and following the prescriptions in Hansen et al., 2014), how the inclusion of that additional source of data uncertainty into a probabilistic inversion via a more complex (non-diagonal) approximation of C d can prevent misinterpretations and artifacts even in the case of environments where the 1D approximation seems reasonable (so, even in the case of settings where the assumption of the absence of significant 3D effects contaminating neighboring 1D soundings seems acceptable).This is something to be investigated in future research; in the present manuscript, we simply evaluate the probability associated with a certain model by using uniquely Equation 2.Moreover, as it will be discussed below, each resistivity component of m will be a function of a specific lithology, allowing the assessment of the probability of occurrence of a lithological category rather than the mere prediction of a resistivity value characterizing a discretizing layer.

Earth and Space Science
10.1029/2023EA003457 ZARU ET AL.

Generation of the Prior Models
The prior distribution can be defined via its realizations.In the case under investigation, the hypothesis is that the subsurface consists of a sequence of two lithologies: peat and clay.It is true that, in some cases, peat might overlay more resistive bedrock (Boon et al., 2008;Comas et al., 2015;Kowalczyk et al., 2017), but our assumption is that, anyhow, even in the case of a shallow bedrock, a clay unit is present due to the rock weathering or alluvial/glacial origins (Osman, 2018).In addition, we assume that, for each sounding location, the associated model can be composed of three lithological units, with no restrictions concerning their vertical sequence.So, for example, a clay unit comprised between a shallower and a deeper peat layer can, in principle, be allowed.Consequently, each 1D model of the prior distribution is discretized with 200 layers with a thickness of 0.1 m, and two interface depths are picked from a uniform distribution with values ranging from 0 to 20 m.This procedure defines the three units.Subsequently, a lithology (either peat or clay) is randomly associated to each of the layers.Therefore, prior sample models-consisting of 200 layers (gathered in 3 lithological units)-are generated.For each layer, it is now possible to match a resistivity value drawn from two distributions describing the ranges of resistivity values expected for each lithology (bar plots in Figure 1a).In order to ensure some additional vertical spatial coherence, the original 200-layer resistivity models are smoothed by means of a 5-layer sliding average filter.In this way, the assumption of independence between neighboring layers is loosened up and a vertical spatial correlation between the 1D model parameters is introduced.The rationale behind this choice is that renouncing the high entropy choice generates more geological plausible 1D resistivity profiles.After this model parameters conditioning, the associated resistivity distributions for the peat and clay units are different from the original ones (solid lines in Figure 1).The 5-layer width for the filter has been chosen by inspecting the effect on the final resistivity vertical profiles: a shorter span would have not satisfactorily reduced the erratic behavior of the original resistivity (thick lines in Figure 1), whereas a filter window larger than 0.5 m would have smeared out the sharp interfaces between the lithologies.
The resistivity ranges in Figure 1 have been deduced from the examples in the literature.Peat soil has a wide resistivity range spanning from tens to a few hundred Ωm, depending on the moisture content, salt concentration, degree of decomposition, and mineral composition (Mohamad et al., 2021;Silvestri, Christensen, et al., 2019;Silvestri, Knight, et al., 2019;Valois et al., 2021).Also the resistivity of clay sediments may vary significantly depending on the soil conditions (in particular, the water content and its salinity) (Cheng et al., 2019;Samouëlian et al., 2005;Solberg et al., 2014).However, for the fresh-water-saturated clay expected in the Alpine peatland, low resistivity values of the order of a few tens Ωm seem to be reasonable (Binley & Slater, 2020).
Therefore, based on these assumptions concerning the possible lithological stratification and the resistivities associated with such geological units, the samples of the prior distribution can be easily generated.Figure 2 displays two examples of the prior distribution models.It might be worth noticing that three lithological units provide additional flexibility since the case of two effective units (e.g., peat covering clay-Figure 2a) is still a feasible option, but, at the same time, it is also possible to experience a clay substrate with a peat intrusion (Figure 2b).Clearly, a very informative prior (e.g., only two units) might be misleading in case of wrong initial assumptions, whereas, to a certain extent, allowing a higher variability in the prior assumptions cannot be detrimental (the only drawback is connected to the higher computational time).
Following this strategy, we populated (and defined) the prior ensemble with 10 5 samples.For each of them, the associated FDEM response has been calculated by taking into account the geometry parameters, frequencies, and elevation of the employed sensor.In this respect, the measured data set consists of 153,621 soundings collected by using the following 11 frequencies: 1, 025,1,525,2,875,5,825,7,775,12,775,15,325,25,525,36,225,63,025,80,225 Hz, with the recording device assumed at 1.0 m height from the surface.During the survey, we used a GEM-2 (Haoping & Won, 2003) in the horizontal co-planar configuration.This device has a receiver-transmitter inter-coil distance of 1.6 m.Performing the standard FDEM processing procedure preliminary to the inversion, we laterally stacked and averaged between adjacent soundings in order to remove the outliers and assess the noise in the data (cf., e.g., Zaru et al., 2023).
For every sounding, the response F(m) of each of the prior realizations is evaluated against the observed data d obs through Equation 2. The posterior distribution is then defined by replicas of the models based on their probability.Hence, the probabilistic inversion consists of an independent extended Metropolis algorithm, which is analogous to the standard extended Metropolis algorithm (Hansen et al., 2016;Mosegaard & Tarantola, 1995) with an infinitely long step-length; this is equivalent to state that every new model proposal from the prior distribution is a completely independent realization: indeed, in our specific case, it is precalculated.The proposed inversion scheme is implemented in the SIPPI toolbox (Hansen et al., 2013a(Hansen et al., , 2013b)).

Results
Figure 3 shows the horizontal slices at different depths (namely: 0.5, 3.0, 5.0, 8.0 m) of: (a) the realization of the posterior distribution with the maximum likelihood (first column on the left); (b) the mean model obtained by averaging the realizations of the posterior models (second column from the left); (c) the probability of encountering peat; and (d) the probability of clay occurrence.Clearly, since we are simply assuming two lithological categories (peat and clay), the two probabilities are complementary: the maximum probability of having peat corresponds to the minimum probability of clay.It is fundamental to clarify that across the present paper, for the sake of conciseness, we use the term maximum likelihood model to indicate the model of the prior ensemble which is characterized by the highest likelihood; strictly speaking, that is not the actual maximum likelihood model since it relates to the prior distribution.However, in this context, and after this caveat, there should be no room left for ambiguities.
From Figure 3, it is evident that peat is mainly present in the northern-central portion of the investigated area and its evolution can be followed from the surface to significant depths (more than 8.0 m).In addition, a shallower peat deposit can be detected in the southern part of the study area.It is worth highlighting the spatial coherence of the result (including the maximum likelihood one) despite each sounding has been inverted separately without introducing any horizontal constraint (again, prior information is provided uniquely via the prior distribution).
In order to verify our result, we compare the probabilistic outcome against the available boreholes and electrical resistivity models obtained by a deterministic inversion.In this respect, Figures 4c-4f shows a vertical section extracted from the (pseudo-)3D volume in Figure 3.All the boreholes were drilled aiming at reaching the clay bottom unit (more difficult to penetrate); however, that was not always possible.So, the depth of the boreholes is supposed to be shallower or match the top of the clay substrate.For sake of clearness, the names of the boreholes reaching the top of the clay unit contain an asterisk (e.g., B26*).From Figures 4c-4f, it is evident that the geophysical reconstruction can detect quite precisely the morphology of the peat unit as it appears from the boreholes.In particular, that is true for the boreholes B26* and B11, in the middle of the section.On the other hand, a possible explanation for the mismatch with the borehole B31* is that it is actually quite off the profile.In fact, the borehole B31* is approximately 4 m apart from the profile.It is worth mentioning that the borehole B26* is the most reliable well as, differently from the others, its core has been extracted and analyzed (Poto et al., 2013).Concerning the much shallower B33, that borehole simply confirms that, at that location, the peat substrate is thicker than 2.5 m, and this evidence is compatible with what is inferred from geophysics.On the contrary, accordingly to the borehole B6*, at that location, the top of the clay unit should be found at 1 m depth, whereas the probable peat thickness (Figure 4e) is supposed to be around 4 m.A plausible justification for such a mismatch will require further investigations; however, because of the noticeable lateral consistency of the geophysical solution, despite the fact that each individual sounding is inverted separately, it is hard to believe in a systematic artifact in the inversion process.In addition, as it will be discussed in the following, similar conclusions about the B6* location can be drawn also by looking at the alternative deterministic reconstruction (Figure 4b).
In general, as it should be, many of the details visible in the maximum likelihood model (Figure 4c) disappear when we check the mean model (Figure 4d) or the peat probability (Figure 4e).This is particularly evident at depth, below the conductive clay, where the model parameters do not have a significant sensitivity with respect to the (noisy) measurements; in other terms, those layers are below the DOI (Depth Of Investigation-e.g., Vignoli et al., 2021;Zhdanov, 2002) and the details retrieved in some specific realizations of the posterior distribution are, indeed, artifacts and are characterized by a low occurrence probability.When we check the data fitting of the different models in terms of

Earth and Space Science
, in which N is the number of data points per sounding (so, in this specific case, N = 11, with 11 being the number of used frequencies), we can see that, not surprisingly, the maximum likelihood model (Figure 4a-red "+" markers) is characterized by responses closer to the observations than the mean model (green "×" markers).Besides, the maximum likelihood model generally fits the data below the optimal threshold χ 2 = 1.Interestingly, also the mean model, which does not have any physical meaning, sensu stricto, is associated to reasonable χ 2 values.
The maximum likelihood solution is equivalent to an inversion obtained by minimizing the data misfit in which the well-posedness of the inversion is enforced by selecting the solution among the prior models (Zaru et al., 2023) rather than via the regularization term.The impact of the prior information is evident as the maximum likelihood reconstruction is characterized by a significant lateral coherence (i.e., the solution turns out to be stable with respect to perturbations in the data-i.e., moving from one sounding to the adjacent one) without an explicit regularization term, and this while the data misfit is comparable with alternative, "more standard," inversions.In this perspective, it might be worth it to compare our probabilistic outcome with the deterministic Occam's inversion (Figure 4b).In doing that, we used the extremely valuable open-source McLachlan et al. (2021b)'s implementation: EMagPy, with a parameterization similar to the one used for the probabilistic inversion.Hence, for example, 30 layers discretize each 1D model; 30 layers seemed to be a good compromise between using a toocoarse discretization (with the risk of introducing an additional, unwanted, regularization parameter: the number of layers) and the practical impossibility of running 200-layer inversions in a reasonable amount of time.In the deterministic approaches, the calculation of the Jacobian is generally the computational bottleneck of the minimization process (Christiansen et al., 2016); even if, in the present research, we have not profusely investigated this aspect, most likely, within the conjugate gradient minimization utilized, the Jacobian's update makes the cost of the deterministic inversion too high to be able to run it for the entire survey (even with a relatively small number of layers).For this reason, we performed the deterministic inversion uniquely along the section highlighted in Figure 3 and shown in Figure 4b.The starting model for each of the 1D deterministic inversions along the investigated profile is the same, with a decreasing resistivity with depth, ranging from around 500 Ωm at the surface to approximately 100 Ωm at the bottom.The regularization weight, λ (Equation 4), is set for all soundings equal to 0.07.This uniform value of λ along the entire profile guaranties an optimal mean data fitting (χ 2 = 1) across the section (cf. the horizontal blue dotted line in Figure 4a), despite, in principle, it should be adjusted in order to ensure that similar condition being fulfilled sounding-by-sounding and not merely on average.Figure 4b shows the deterministic reconstruction, whereas the associated data fitting levels can be checked in Figure 4a.In a few sounding locations, the deterministic inversion has difficulties in fitting the data to a reasonable level.
Possible reasons for that could be related: to local minimization issues of the objective functional, and/or to the inconsistency between the prior information enforced by the stabilizer, which promotes smoothly varying solutions, not really compatible with the geology characterized by quite a sharp peat-clay interface (sometimes occurring at shallow depths).In this perspective, a sounding-by-sounding optimization of the λ value might have been useful even if extremely time and computationally consuming.The deterministic inversion has a data misfit that is largely overlapping with the one of maximum likelihood model; still, it is more laterally erratic, while, not surprisingly, favors the retrieval of vertically elongated structures.On this matter, the comparison, for example, of the portion of the sections between B26* and B11 as inferred by the two alternative inversion schemes, demonstrates the impact of the different prior information incorporated in the probabilistic approach and in the deterministic regularization: the complex information formalized via the prior distribution samples allows the retrieval of a spatially consistent morphology of the peat unit, whereas this is prevented by the smooth regularization (even in case of similar data fitting).Earth and Space Science 10.1029/2023EA003457 In addition, the good correspondence between the probabilistic result and the deterministic reconstruction in the vicinity of borehole B6* questions the reliability of the direct measurement of the peat thickness deduced by the drilling.
In terms of inversion performances, the probabilistic inversion basically requires only the forward calculations of the samples of the prior distribution, which can be performed independently and, so, are massively and easily parallelizable.The calculation of the conditional probability of Equation 2 is, then, straightforward and, in the specific case of the entire Danta's data set (consisting of around 150,000 soundings-a size comparable with an airborne survey; cf., e.g., Christiansen et al., 2016;Viezzoli et al., 2009), with a prior ensemble of 100,000 samples, took approximately 1 hr (on a quite standard laptop equipped with an Intel i7-10750H/2.60GHz).Clearly, on top of that, for a fair comparison, the time used for generating the prior samples and the associated measurements must be considered (for the 100,000 samples, on the same laptop, it took around one additional hour).However, it should be stressed that for similar geological settings (and the same configurations of the geophysical instrumentation), the prior models do not need to be recalculated every time and can be generated once and for all.On the other hand, at least with the parameterization we used, the deterministic inversion of the entire Danta's data set was not feasible (it would have needed weeks of calculation), and, for each of the several 2D profiles it was tested on (like in Figure 4), required more than one day.Furthermore, it might be relevant to stress that even if populating the prior distribution might be difficult, similar troubles might be encountered with deterministic inversion schemes, not only in the definition of the most appropriate stabilizer, but also in the selection of the optimal regularization weight and in addressing all the issues connected with the optimization of the objective functional (including, e.g., the definition of a suitable starting model).It is worth stressing once more that the high performances characterizing the proposed probabilistic approach relate to the proposal strategy: since all the evaluated models are mutually independent, they can be precalculated in advance (in parallel) (Hansen, 2021a).

Discussion
Even if, in presenting the inversion results we referred repeatedly to the mean map (Figure 4d), clearly, the actual probabilistic outcome (the posterior PDF) contains much more information than the mean model and its associated standard deviation.In this respect, it is important to stress, that, also because of the non-linearity of the forward modeling, the resulting PDF might be non-Gaussian.Therefore, the possible multi-modality would make the reduction of the probabilistic result to its mean and standard deviation misleading.In this respect, four sounding locations (Figure 5a) have been selected to show the risk of taking just the mean as an immediately interpretable result; in particular, in Figures 5b-5e, the posterior resistivity PDFs associated with those four locations are plotted and compared against the mean resistivity vertical profile (red solid line), together with the corresponding standard deviation band (between the green solid lines).Consistently with the previous caveat, the resistivity PDFs present, at depth (specifically, where the occurrence of clay is more probable), a bi-modal pattern, which, in turn, makes the mean resistivity move between the probability maxima.
Hence, panels in Figure 5 confirm once more the importance of a probabilistic approach as, in this multi-modal situation, also deterministic approaches (e.g., orange dashed line in Figure 5b) could merely provide a local assessment of the uncertainty (based on the Jacobian) associated with the unique retrieved solution.And, of course, the same risk may take place if we draw our conclusions merely based on the mean model.Also, it is worth noticing how the probability of peat occurrence (magenta solid line) drops rapidly when the gyttja starts to be encountered (on the left side of Figure 5b, the detailed lithological description of borehole B26* can be found).probability density functions (PDFs) as a function of depth for each of the four selected locations (respectively associated with soundings no.26805, 50,000, 100,000, and 150,000); for each of the selected soundings, the posterior PDF is compared against the corresponding mean model resistivity (solid red line) and the relative standard deviation limits (the band between the solid green lines), and, for completeness, also with the maximum likelihood model (black solid line) and the peat probability (solid magenta line, whose x-axis is shown on the top of each panel).For the sounding 26,805, the deterministic inversion is also available and, consistently, is plotted in panel (b) as a dashed orange line.Sounding 26,805 was selected since it is the closest to the borehole B26* whose stratigraphy is known and appears alongside panel (b).
In Figure 5, the vertical resistivity profile corresponding to the maximum likelihood solution (black solid line) is also plotted for each of the picked sounding locations.The maximum likelihood solutions are those (within the prior ensemble) better fitting the data, not the most probable.Indeed, the maximum likelihood model is the maximizer of the conditional PDF, p(d obs |m), in Equation 2, and not of the posterior PDF, p(m|d obs ).In the presented framework, since the prior is defined uniquely via its samples and there is no closed form for it, the calculation of the most probable model is not possible.
Despite the maximum likelihood model being simply the model better fitting the observations, in the present scheme, to some extent, it is still a regularized solution as the prescriptions infused by the prior distribution impose a non-arbitrary final model.As mentioned before, this is clear, for example, from Figure 4c: even though neighboring soundings (characterized by similar measurements) are independently inverted, the corresponding results are similar (demonstrating a continuity of the models with respect to the data).And this is particularly remarkable because no lateral constraints have been enforced.In other circumstances, it might be useful to constrain the inversion also horizontally.An extension of the present scheme capable of considering the expected lateral continuity might consist of adding a term taking into account the probability of a 1D model being different from some reference model (possibly, the reference model can result from the previously considered adjacent 1D sounding inversion).Similar attempts are described, for example, in Ardid et al. (2021) and Zaru et al. (2023).

Conclusions
We presented a workflow for the probabilistic petrophysical 1D inversion of FDEM measurements and tested it on a data set consisting of 11 frequencies collected throughout an Alpine peatland in northern Italy.The retrieved posterior distribution has been analyzed in terms of maximum likelihood and mean models, and via the associated peat (and clay) occurrence probability.It turned out that the peat probability can effectively reconstruct the geometry of the peat layer and satisfactorily match the additional ancillary pieces of evidence about the investigated site.In particular, the probabilistic results have been compared against the available boreholes and their capability was confirmed at least for the deepest (and most reliable) direct investigations.Furthermore, the probabilistic scheme was also verified against the most standard deterministic Occam's inversion implemented in a very well-known open-source software package.The deterministic and the probabilistic inversion show remarkable differences due to the different prior information formalized in the two schemes (smooth vertical constraints in the deterministic inversion; three lithological units with specific resistivity ranges in the probabilistic approach).We also highlight that the adopted probabilistic approach can be extremely efficient paving the way to quasi-real-time inversion in the field and on-the-fly optimization of the survey plans.
Future developments should involve the assessment of the modeling error and the inclusion of such an inevitable source of uncertainty into the inversion scheme as has already been done for the time-domain counterpart (Bai et al., 2021).

Figure 1 .
Figure 1.Resistivity distributions associated with each lithology.Panel (a) displays, as green bars, the probability density distribution (PDF) of having a peat layer with a specific resistivity, and, as yellow bars, the corresponding PDF for the clay layers.The solid lines show the corresponding distributions after a 5-layer vertical moving average filtering of the samples of the prior.In panel (b), the corresponding cumulative density functions (CDFs) are shown: the solid lines are related to the CDF obtained after the 5-layer vertical filtering and highlight that resistivity values lower than 200 Ωm are almost certainly associated with clay, whereas layers with resistivities higher than 400 Ωm, most likely, are peat.Resistivities around 270 Ωm have the same probability of being clay or peat.

Figure 2 .
Figure 2. Samples of the prior distribution.Each prior model consists of 200 layers characterized by different resistivity values (thin green or orange solid lines) gathered in three lithological units (of either clay or peat).The original 200 resistivities have been subsequently filtered to obtain 1D model samples with some vertical correlation (thick blue or red solid lines).The mean resistivity calculated from the filtered values over each lithological unit is shown as a black thin straight segment.Panel (a) displays a sample with the following lithological sequence: peat-clay-clay; consequently, the mean resistivities for the two deepest units are almost indistinguishable (their interface occurs at around 5.5 m depth).Panel (b) illustrates another example of a prior sample characterized by a peat unit sandwiched between two clay units.

Figure 3 .
Figure 3. Horizontal slices of the inversion results of the frequency-domain electromagnetic (FDEM) data collected across the Danta peatland (Italy).Each row displays, from left to right: (i) the resistivity solution with the maximum likelihood; (ii) the mean model of the samples of the posterior probability; (iii) the probability of peat; and (iv) the probability of clay.Rows (a)-(d) show the horizontal slices of those models and probabilities, respectively, at depth 0.5, 3.0, 5.0, and 8.0 m.In the column 3, on top of the peat probability, the location of the profile in Figure 4 is shown (in magenta), together with the position of the related boreholes (green circles).

Figure 4 .
Figure 4. Inversion results.Panel (b) shows the fully nonlinear deterministic Occam's inversion performed with the software package EMagPy (McLachlan et al., 2021b).Panels (c)-(f) display the results of the probabilistic petrophysical inversion; in particular, panel (c) shows the maximum likelihood resistivity solution, whereas panel (d) exhibits the resistivity section obtained by averaging the resistivity samples of the posterior distribution.Panels (e) and (f) show respectively the probability of peat and clay occurrences for each of the 200 layers discretizing each 1D model.Panel (a) shows the χ 2 value obtained by comparing, against the observed data, the frequency-domain electromagnetic (FDEM) responses of: (i) the deterministic solution (blue dots); (ii) the maximum likelihood model (red crosses); and (iii) the mean model (green "×" markers); the horizontal dotted lines are simply the mean χ 2 value across the section for the three models in panels (b)-(d).The parameterization for the two kinds of inversion is different: 200 layers are used for the probabilistic inversion; 30 for the deterministic one.

Figure 5 .
Figure 5.Comparison of individual soundings.Panel (a) shows the locations of the selected frequency-domain electromagnetic (FDEM) soundings and the position of the borehole B26*.Panels (b)-(e) display the electrical resistivity probability density functions (PDFs) as a function of depth for each of the four selected locations (respectively associated with soundings no.26805, 50,000, 100,000, and 150,000); for each of the selected soundings, the posterior PDF is compared against the corresponding mean model resistivity (solid red line) and the relative standard deviation limits (the band between the solid green lines), and, for completeness, also with the maximum likelihood model (black solid line) and the peat probability (solid magenta line, whose x-axis is shown on the top of each panel).For the sounding 26,805, the deterministic inversion is also available and, consistently, is plotted in panel (b) as a dashed orange line.Sounding 26,805 was selected since it is the closest to the borehole B26* whose stratigraphy is known and appears alongside panel (b).