Coupled land surface–subsurface hydrogeophysical inverse modeling to estimate soil organic carbon content and explore associated hydrological and thermal dynamics in the Arctic tundra

. Quantitative characterization of soil organic carbon (OC) content is essential due to its signiﬁcant impacts on surface–subsurface hydrological–thermal processes and microbial decomposition of OC, which both in turn are important for predicting carbon–climate feedbacks. While such quantiﬁcation is particularly important in the vulnera-ble organic-rich Arctic region, it is challenging to achieve due to the general limitations of conventional core sampling and analysis methods, and to the extremely dynamic na-ture of hydrological–thermal processes associated with an-nual freeze–thaw events. In this study, we develop and test an inversion scheme that can ﬂexibly use single or multiple datasets – including soil liquid water content, temperature and electrical resistivity tomography (ERT) data – to estimate the vertical distribution of OC content. Our approach relies on the fact that OC content strongly inﬂuences soil hydrological–thermal parameters and, therefore, indirectly controls the spatiotemporal dynamics of soil liquid water content, temperature and their correlated electrical resistivity. We employ the Community Land Model to simulate non-isothermal surface–subsurface hydrological dynamics from the bedrock to the top of canopy, with consideration of land surface processes (e.g., solar radiation balance, evapotranspi-ration, snow accumulation and melting) and ice–liquid water phase transitions. For inversion, we combine a deterministic and an adaptive Markov chain Monte Carlo (MCMC) optimization algorithm to estimate a posteriori distributions of desired model parameters. For hydrological–thermal-to-geophysical variable transformation, the simulated subsurface temperature, liquid water content and ice content are explicitly linked to soil electrical resistivity via petrophysical and geophysical models. We validate the developed scheme using different numerical experiments and evaluate the inﬂuence of measurement errors and beneﬁt of joint inversion on the estimation of OC and other parameters. We also quantify the propagation of uncertainty from the estimated parameters to prediction of hydrological–thermal responses. We ﬁnd that, compared to inversion of single dataset (temperature, liquid water content or apparent resistivity), joint inversion of these datasets signiﬁcantly reduces parameter uncertainty. We ﬁnd that the joint inversion approach is able to estimate OC and sand content within the shallow active layer (top 0.3 m of soil) with high reliability. Due to the small variations of temperature and moisture within the shallow permafrost (here at about 0.6 m depth), the approach is unable to estimate OC with conﬁdence. However, if the soil porosity is functionally related to the OC and mineral content, which is often observed in organic-rich Arctic soil, the uncertainty of OC estimate at this depth remarkably decreases. Our study documents the value of the new surface–subsurface, deterministic–stochastic inversion approach, as well as the beneﬁt of including multiple types of data to estimate OC and associated hydrological–thermal dynamics.

Abstract.Quantitative characterization of soil organic carbon (OC) content is essential due to its significant impacts on surface-subsurface hydrological-thermal processes and microbial decomposition of OC, which both in turn are important for predicting carbon-climate feedbacks.While such quantification is particularly important in the vulnerable organic-rich Arctic region, it is challenging to achieve due to the general limitations of conventional core sampling and analysis methods, and to the extremely dynamic nature of hydrological-thermal processes associated with annual freeze-thaw events.In this study, we develop and test an inversion scheme that can flexibly use single or multiple datasets -including soil liquid water content, temperature and electrical resistivity tomography (ERT) data -to estimate the vertical distribution of OC content.Our approach relies on the fact that OC content strongly influences soil hydrological-thermal parameters and, therefore, indirectly controls the spatiotemporal dynamics of soil liquid water content, temperature and their correlated electrical resistivity.We employ the Community Land Model to simulate nonisothermal surface-subsurface hydrological dynamics from the bedrock to the top of canopy, with consideration of land surface processes (e.g., solar radiation balance, evapotranspiration, snow accumulation and melting) and ice-liquid water phase transitions.For inversion, we combine a deterministic and an adaptive Markov chain Monte Carlo (MCMC) optimization algorithm to estimate a posteriori distributions of desired model parameters.For hydrological-thermal-togeophysical variable transformation, the simulated subsurface temperature, liquid water content and ice content are explicitly linked to soil electrical resistivity via petrophysical and geophysical models.We validate the developed scheme using different numerical experiments and evaluate the influence of measurement errors and benefit of joint inversion on the estimation of OC and other parameters.We also quantify the propagation of uncertainty from the estimated parameters to prediction of hydrological-thermal responses.We find that, compared to inversion of single dataset (temperature, liquid water content or apparent resistivity), joint inversion of these datasets significantly reduces parameter uncertainty.We find that the joint inversion approach is able to estimate OC and sand content within the shallow active layer (top 0.3 m of soil) with high reliability.Due to the small variations of temperature and moisture within the shallow permafrost (here at about 0.6 m depth), the approach is unable to estimate OC with confidence.However, if the soil porosity is functionally related to the OC and mineral content, which is often observed in organic-rich Arctic soil, the uncertainty of OC estimate at this depth remarkably decreases.Our study documents the value of the new surface-subsurface, deterministic-stochastic inversion approach, as well as the benefit of including multiple types of data to estimate OC and associated hydrological-thermal dynamics.

Introduction
Soil organic carbon (OC) and its influence on terrestrial ecosystem feedbacks to global warming in permafrost regions are particularly important for the calculation of global carbon budget and prediction of future climate variation.
Warmer air temperature leads to permafrost degradation, which is expected to enhance decomposition of huge pools of previously frozen OC, releasing carbon dioxide and methane to the atmosphere, and enhancing global warming (Koven et al., 2011;Schaphoff et al., 2013;Schuur et al., 2015).In that context, accurate estimation of OC content stored in both the active layer and permafrost is crucial for investigation of carbon stocks that expose microbial decomposition.
Predictive understanding of ecosystem feedbacks to climate in permafrost regions requires quantitative knowledge of surface-subsurface hydrological-thermal dynamics, which in turn are strongly governed by the hydrologicalthermal properties of soil OC (Jafarov and Schaefer, 2016).In particular, there are dramatic differences between thermal and hydraulic properties of OC and mineral soil, both of which typically co-exist in shallow permafrost systems.OC's thermal conductivity (e.g., λ OC, dry = 0.05 W mK −1 ) is significantly lower than that of mineral soil (e.g., λ sand = 8.4 W mK −1 ) (Farouki, 1981).By contrast, its heat capacity is higher than that of mineral soil.Considering hydrological properties, the hydraulic conductivity of OC is higher and its capillary pressure is smaller than mineral soil (Lawrence and Slater, 2008).In addition, while mineral soil porosity typically ranges from 0.4 to 0.6, the porosity of OC soil is usually greater than 0.8.Due to its low thermal conductivity, a top OC layer can behave as an insulator that reduces the magnitude of heat and energy exchange between the atmosphere and deeper soil (e.g., Hinzman et al., 1991;Rinke et al., 2008).Nicolsky et al. (2007) and Jafarov and Schaefer (2016) reported that inclusion of vertical OC content profile into a land surface model can considerably improve prediction of subsurface moisture, temperature and carbon dynamics.However, our ability to measure or estimate the distribution of OC is currently challenging, which inhibits accurate model prediction.
OC content is usually measured from core samples, which are collected from field sites and then analyzed in the laboratory (e.g., Kern, 1994).While this method is relatively accurate, it is labor intensive and typically limited in spatial coverage.Because OC and mineral content largely influence hydrological-thermal parameters (i.e., thermal conductivity, heat capacity, hydraulic conductivity and retention curve; see Appendix A), they are the main soil properties that control the subsurface hydrological-thermal dynamics.As a result, OC and mineral content can be potentially obtained by inverting observations of hydrological-thermal state variables (i.e., soil liquid water content and soil temperature) and their correlated observables (e.g., electrical resistivity).However, so far there has been no effort using this approach to indirectly estimate these soil properties.
Geophysical methods hold potential for characterizing the subsurface in permafrost regions as well as their associated physical, hydrological and thermal processes.Geophysical techniques offer an advantage over conventional point measurement techniques because they provide spatially exten-sive information in a minimally invasive manner (e.g., Hubbard and Rubin, 2005).For example, Arcone et al. (1998) and Chen et al. (2016) used ground-penetrating radar (GPR) to characterize the depth of the permafrost table.Hinkel et al. (2001) used GPR to estimate thaw depth, to recognize ice wedges and ice lenses, and to locate the organic-mineral soil interface.Schwamborn et al. (2002) combined seismic and GPR data to investigate the stratigraphy of both frozen and unfrozen parts of Lake Nikolay.Lewkowicz et al. (2011) and You et al. (2013) employed electrical resistivity tomography (ERT), ground temperature monitoring, frost table probing and coring to detect the permafrost depth.Hauck et al. (2011) developed a four-phase model of soil matrix, ice, liquid and air and used it to estimate soil liquid and ice content from combined ERT and seismic measurements in the Swiss Alps.Hubbard et al. (2013) combined lidar data with multiple geophysical (ERT, GPR, electromagnetic) and point measurements to characterize active-layer thickness and permafrost variability in a large area.
In spite of the potential benefits offered by geophysical data for characterizing permafrost systems, geophysical inversion approaches typically suffer from several challenges.First, inversion methods are often ill-posed due to the fact that geophysical observables are sensitive to different soil properties.Secondly, inversion approaches often typically require petrophysical models to link the geophysical observables with the property of interest.Finally, there are differences between the geophysical support scale and the scale of the imaging target (Hubbard and Linde, 2011).In order to take advantage of information inherent in geophysical signatures and minimize the non-uniqueness challenges described above, many recent studies have explored the value of coupled hydrogeophysical inversion frameworks for estimating soil properties (e.g., Johnson et al., 2009;Huisman et al., 2010;Irving and Singha, 2010;Kowalsky et al., 2011;Pollock and Cirpka, 2012;Busch et al., 2013;Herckenrath et al., 2013;Camporese et al., 2015;Tran et al., 2013Tran et al., , 2016)).In these studies, the hydrological and geophysical models are coupled together so that geophysical data are used to estimate soil properties that control the subsurface hydrologicalthermal dynamics.Of the geophysical techniques commonly used for monitoring the shallow subsurface, ERT is increasingly common because it can autonomously provide 2-or 3-D time-lapse measurements with a relatively high spatial resolution, is sensitive to properties influencing hydrologicalthermal dynamics and is particularly suitable for field deployment over a long period of time.As a result, we use ERT data in this study.
Most coupled hydrogeophysical inversion approaches developed to date are not adequate for investigating permafrost systems due to several gaps.Developed methods have only been applied to terrestrial systems without consideration of the significant dynamics associated with the freeze-thaw transition.Developed coupled hydrogeophysical inversion approaches have also not yet incorporated The Cryosphere, 11, 2089Cryosphere, 11, -2109Cryosphere, 11, , 2017 www.the-cryosphere.net/11/2089/2017/surface-subsurface interactions (e.g., evapotranspiration, energy balance, plant water uptake).Finally, while a few studies have used soil-vegetation-atmospheric transfer (SVAT) models to qualitatively interpret geophysical data (e.g., Mc-Clymont et al., 2013), to date no study has coupled SVAT and geophysical models and data to improve property estimation.
Building on recent advances in the use of electrical methods in the permafrost (e.g., Minsley et al., 2015;Dafflon et al., 2017) as well as coupled hydrogeophysical inversion approaches described above, this study focuses on the development of an inverse approach that uses single or multiple datasets (soil liquid water content, soil temperature and electrical resistivity) to estimate OC content, which is a main factor that governs the subsurface hydrological-thermal dynamics.Our approach advances and couples several algorithms.We use a SVAT model known as the Community Land Model (CLM4.5;Oleson et al., 2013) to simulate water, heat and energy exchange from the bedrock to the top of the canopy.The model considers most of the land surface processes, iceliquid phase change and surface-subsurface hydrologicalthermal dynamics.For parameter estimation, we combined deterministic and stochastic optimization algorithms to concurrently obtain the best parameter estimates and their associated uncertainties.The deterministic optimization algorithm is employed to estimate the initial parameter set and covariance matrix of the proposal distribution.For the stochastic optimization, we used an advanced Markov chain Monte Carlo (MCMC) method known as delayed rejection adaptive Metropolis (DRAM; Haario et al., 2006).With this implementation of this adaptive MCMC algorithm, we expect to obtain the a posteriori probability density function (PDFs) of the desired model parameters more quickly than with the traditional MCMC technique.For hydrological-thermal-togeophysical transformation, we explicitly consider the dependence of the soil electrical resistivity on the soil iceliquid water content and soil temperature via petrophysical and forward geophysical models.
This study advances capabilities to estimate and understand the controls of OC on hydrological and thermal properties through developing a hydrological-thermal-geophysical inversion scheme and through exploring its potential to estimate the vertical distribution of OC and mineral content at several depths within a representative synthetic Arctic soil column.Herein, we use synthetic studies to (1) evaluate the relationship between the measurement error and uncertainties of parameter estimates; (2) examine the improvement in parameter estimation offered by including various datasets in the inversion, including apparent resistivity data; (3) investigate how OC estimation changes if the mineral and petrophysical parameters are unknown; (4) explore how parameter estimation changes when soil porosity functionally correlates with the OC and mineral content; and (5) investigate the uncertainty propagation from the OC and mineral content to the hydrological-thermal prediction.
The paper is organized as follows.Section 2 describes the development of the hydrological-thermal-geophysical inversion scheme.Section 3 analyzes and discusses the results of different synthetic experiments.Summary and concluding remarks are provided in Sect. 4.

Methodology
Generally, the joint hydrological-thermal-geophysical inversion scheme developed in this study (Fig. 1) includes two main components: (1) a forward coupled hydrologicalthermal-geophysical model that generates the subsurface state variables (i.e., ice-liquid water content and temperature) and then uses these variables to infer the apparent resistivity using a set of petrophysical formulas and a forward electrical resistivity model (Fig. 1a), and (2) a combined deterministic-stochastic optimization algorithm to estimate the PDFs of desired model parameters (p) -which include the soil OC content vertical profile (scenarios 1 to 9), sand content vertical profile (scenarios 8 and 9) and petrophysical parameters (scenarios 8 and 9) (see Table 2) -by minimizing the misfit between measured and simulated data.It is worth noting that the scheme is developed so that single (e.g., soil temperature, liquid water content or apparent resistivity) or multiple datasets can be used for inversion.

Hydrological-thermal model
In this study, we employed the CLM4.5 model (hereafter referred to as "CLM"), which can effectively simulate different land surface energy balance and surface-subsurface hydrological-thermal processes (Oleson et al., 2013).CLM represents horizontal heterogeneity using multiple parallel soil and snow columns having different land use and plant function types.The lateral flow between the soil columns is not accounted for in CLM.The model simulates the freezethaw dynamics by considering two phases of water: liquid and ice.The rate of phase change depends on the energy excess (for the ice-to-liquid transition) or deficit (for the liquidto-ice transition) from the soil temperature to the freezing temperature.Given CLM's ability to simulate different hydrological-thermal processes in cold regions, we found it suitable for Arctic soil column simulations.The minimum requirements for the top boundary conditions in CLM include precipitation, incident solar, air temperature and wind speed.The land use and plant type information can be provided by users or extracted from the available model database.
CLM assumes that soil is a mixture of three soil types, namely, OC, sand and clay.It calculates the soil hydrological-thermal parameters based on the content (fraction) of these soil types and their corresponding hydrological-thermal properties (see Appendix A for more detailed information on these relationships).In CLM, the soil OC content (% OC) is defined as % OC = ρ OC /ρ max OC ,  2).The scheme permits flexibly using single or multiple types of data for inversion.The forward coupled hydrological-thermal-geophysical model (a) is iteratively executed in both deterministic and stochastic inversion stages.
Table 1.Petrophysical parameters and soil properties information used for synthetic simulation.Petrophysical parameters Na and β − Cl are obtained from Minsley et al. (2015).Soil OC and sand content are based on the core sample analysis at the site near Barrow, AK.The sand content is the percentage of sand in the mineral mixture, which is calculated as 100 minus OC content.The soil porosity is independent from soil OC and sand content for scenarios 1-8, while it is calculated from these properties in scenario 9.

Petrophysical parameter
Soil properties 9.6487 × 10 4 OC and sand (%) mixture) content) (%) in which ρ OC is the soil OC density (kg m −3 ) and ρ max OC is the maximum soil OC density (ρ max OC = 130 kg m −3 ), which is the standard bulk density of peat (Oleson et al., 2013).The mineral content is determined as % mineral = 100 − % OC.CLM further assumes that mineral includes only sand and clay.As a result, in the inversion scheme, we only need to estimate the soil OC content and sand content in mineral (hereafter referred to as the sand content).The clay content is obtained by subtracting the sand content from the mineral content.
For more detailed exploration of the vertical variability of subsurface properties and associated hydrological-thermal dynamics, we increased the default number of soil layers in the CLM from 15 to 32 layers and defined the depth of layer i (z i ) as z i = 0.025 e 0.17(i−0.2) − 1 . (1) Of these 32 layers, CLM assumes that the 5 bottom layers are bedrock layers.Hydrological dynamics is simulated only in the top 27 soil layers, while thermal dynamics is simulated in all 32 layers.Equation ( 1) was used to ensure that the layer thicknesses near the soil surface are thinner than those near the bottom (as shown in Fig. 3) in order to capture the important hydrological and thermal dynamics in the topsoil active layers.Moreover, in order to explore how the soil porosity influences the estimation of soil OC and sand content, we modified the CLM to consider two cases: (1) the soil porosity profile was fixed and independent from the soil OC and sand content (see scenarios 1-8 in Table 1), and (2) the soil porosity was calculated from the OC and sand content as the default in the CLM (see scenario 9 in Table 1) as below (Lawrence and Slater, 2008): in which is the soil porosity, and min and OC are the porosity of mineral and OC, respectively.In the CLM, the OC porosity is given as OC = 0.9, and the mineral porosity is calculated from sand fraction as min = 0.489 − 0.00126(%sand). (3) The dependencies of soil thermal conductivity, heat capacity and thermal diffusivity on liquid water saturation, OC and sand content are shown in Fig. 2.This figure was obtained from calculations using equations in Appendix A in which the soil porosity was considered in two cases: (1) fixing at 0.7 (panels a-c) and (2) calculating from the OC and sand content (panels d-f).The figure shows that the variation of soil thermal properties with respect to the OC content, sand content and liquid water saturation is similar for both cases.
When the OC fraction increases from 0 to 100 %, the soil thermal conductivity decreases, and the soil heat capacity slightly increases.By contrast, higher sand fraction leads to higher thermal conductivity and slightly lower heat capacity.These relationships are expected, given that OC has a considerably smaller thermal conductivity and a slightly higher heat capacity than sand.The figure also shows that both soil thermal conductivity and heat capacity significantly increase with increasing liquid water saturation.This is also reasonable, as the thermal conductivity and heat capacity of liquid water are much higher than those of air.The thermal diffusivity is defined as the ratio between the thermal conductivity and heat capacity.The figure indicates that the diffusivity increases when the OC decreases and sand content increases.
Comparing the two cases shows that, when the soil porosity depends on OC and sand content, the soil thermal properties change in larger ranges with the variation of OC content, sand content and liquid water saturation.It is because, while the soil porosity is fixed at 0.7 in the first case (panels a-c), it varies from 0.36 (when soil is 100 % sand) to 0.9 (when soil is 100 % OC) in the second case (panels d-f).Bewww.the-cryosphere.net/11/2089/2017/The Cryosphere, 11, 2089-2109, 2017 Table 2. Information on parameter, "observation" data and inversion results for nine scenarios.For scenarios 1-7, the soil OC content (three parameters) at z = 0.1, 0.3 and 0.6 m was estimated.For scenarios 8 and 9, the OC content at z = 0.1, 0.3 and 0.6 m; the sand content in the mineral mixture at z = 0.3, 0.6 and 1 m; and the petrophysical parameters m and α were estimated (eight parameters).The best-estimated parameters and their uncertainties are represented by the mean and standard deviation of the MCMC samples.cause soil thermal properties strongly depend on soil porosity (see Eqs. A2 and A8 in Appendix A), together with the OC and sand content, the variation of porosity in the second case leads to rapid change of the soil thermal properties and, therefore, of the subsurface hydrological-thermal dynamics.

Petrophysical and geophysical transformation
In our inverse scheme, we link the output of the hydrological-thermal simulation described above (soil iceliquid water saturation and temperature) to soil electrical conductivity using Archie's law (Archie, 1942): in which φ is the porosity; S w is the liquid water saturation in the pore space; m and n are the cementation and saturation indexes, respectively; and σ s is the soil electrical conduction, which was fixed at σ s = 0.005 S m −1 in this study (Table 1).
It is worth noting that the reduction of porosity due to ice content in this study was not considered.How ice content influences Archie's equation will be considered in future research.
The water electrical conductivity (σ w ) is calculated from the concentration of all ions in water as (Minsley et al., 2015) in which β i and z i are the ionic mobility and valence of the ith ion, respectively.Similar to Minsley et al. (2015), we assumed that Na + and Cl − are the two main ions in this synthetic study.F c is Faraday's constant.C i is the concentration of the ith ion, which depends on the ice / liquid water fraction as in which S fi and S fw are, respectively, the fraction of ice and liquid in ice-liquid water (S fi + S fw = 1); α varies from 0 to 1, which is the coefficient accounting for the reduction of soil water salinity when liquid water saturation decreases.
A larger α implies a larger increasing rate of ion concentration with decreasing liquid water fraction.The concentrations of ions in the ice-free water (C i(S fi =0) ) can be obtained from samples in the summer season.α and m are the two most important parameters that control the relationship between geophysical and hydrologicalthermal variables.We estimated them by inverting soil moisture, temperature and geophysical data in scenarios 8 and 9 (see Table 2).The effect of soil temperature (T ) on the soil electrical conductivity is formulated as (Hayley et al., 2007) σ T = σ (0.018 × (T − 25) + 1). (7) The linkage between soil electrical conductivity and the apparent resistivity is established by the electrical forward model.In this study, we used the forward model of the Boundless Electrical Resistivity Tomography (BERT) package, developed by Rücker et al. (2006), which numerically solves Poisson's equation using the finite-element method in a three-dimensional arbitrary topography.For more detailed information on this model, we refer the reader to Rücker et al. (2006).

Stochastic and deterministic parameter estimation
In this section, we present a combination approach of deterministic and stochastic optimization algorithms to estimate the model parameters p and their uncertainties.The stochastic optimization algorithm relies on the Bayesian inference and DRAM MCMC technique.The deterministic optimization algorithm was used to approximate the initial set of model parameters and initial covariance matrix of the proposal distribution for stochastic optimization.Consequently, the estimated parameters are more rapidly obtained than only using a single stochastic algorithm with arbitrary initial parameters.Moreover, the use of the DRAM stochastic optimization algorithm allows us to sequentially update the proposal covariance matrix and perform multiple tries to improve the acceptance rate.This algorithm has proven to be more efficient than the commonly used MCMC Metropolis-Hasting method.

Bayesian inference
In the stochastic parameter estimation, the objective is to find the a posteriori probability distribution P (p|Y ) of parameters p conditioned on the measurements Y from which we can extract the best-estimated parameters and their uncertainties.Based on Bayes' rule, this a posteriori distribution is formulated as follows: in which p (p) is the a priori parameter distribution of parameter p and p (Y |p) is the likelihood function.Assuming that the error residuals are uncorrelated, the likelihood function can be written as where f y i (y i |p) denotes the PDF of measurement y i at time t i given the model parameters p.If we further assume the error residuals (difference between modeling and measurement) to be normally distributed, then f y i (y i |p) can be written as and Eq. ( 9) becomes in which y i (p) is the model response at time t i and σ 2 i is the variance of measurement error at time t i .Intuitively, σ 2 i works as an inverse-weighted factor of the contribution of measurement ŷi to the a posteriori distribution p (Y |p).A measurement with a higher variance of measurement error has a smaller contribution to constructing the parameter a posteriori distribution.In addition, for joint inversion, σ 2 i helps to removes the influence of measurement units of different data types.

Delayed rejection adaptive Metropolis Markov chain Monte Carlo method
Once the a posteriori density distribution p (p|Y ) of the model parameters is defined, we need to determine its statistical properties (e.g., mean, covariance).However, due to the nonlinearity of the dynamic model, it is usually difficult to analytically obtain these properties.In that respect, the Monte Carlo methods can be used to generate samples from this a posteriori distribution and then calculate these properties.We employed the DRAM method that was improved from the Metropolis-Hasting MCMC method for this purpose.Basically, this method is a combination of the adaptive Metropolis and delayed rejection algorithm and briefly presented as follows: Metropolis-Hasting: Given the current parameter set p k at iteration k, the candidate for the next move (p k+1 ) from the current value is generated from a proposal distribution q 1 (p k , p k+1 ).The acceptance ratio is calculated as below: If the second try is rejected, the third try can be generated and so on.The number of tries is specified by users.
Adaptation: One of the key limitations of the MCMC technique is the selection of the proposal distribution model.In adaptive Metropolis, the proposal distribution is assumed to be Gaussian centered at the current sample N (p k C k ) with the covariance matrix C k adapted from the previous samples as In Eq. ( 14), s d is the scaling parameter, which depends on the length (d) of the estimated parameter vector p, which is set to s d = 2.4 2 /d; ε > 0 is a very small constant to inhibit C k from becoming singular; and I d signifies the d-dimensional identity matrix.
The assessment of convergence of the MCMC chain is analyzed by Geweke's criterion, which compares the means and variances of the beginning and end segments of the chain as below: where a denotes the beginning interval, which was selected as the first 10 % of the chain, and where b denotes the end interval, which was selected as the last 50 % of the chain.p i,a and p i,b are, respectively, the mean of the ith parameters of segments a and b; n a and n b are the number of samples in a and b segments; and s i,a and s i,b are their corresponding consistent spectral density estimates at zero frequency.The chain is considered to be converged if the G i score is within the 95 % interval of the standard Gaussian distribution (−1.96 ≤ z i ≤ 1.96).

Deterministic optimization for approximating starting parameters and proposal distribution
The speed of convergence of the MCMC optimization algorithm strongly depends on the initial model parameter p 0 and initial proposal distribution q 1 .In order to reduce the number of iterations needed to obtain the a posteriori distribution of model parameters, we used the local optimization Nelder-Mead simplex method to approximate the starting model parameter and initial covariance matrix of the proposal distribution.The starting point of the DRAM is the best-estimated set of model parameters obtained by the Nelder-Mead method.
The covariance matrix of the proposal distribution is assumed to be similar to that of the model parameters, which are locally calculated at the optimal solution obtained by the Nelder-Mead method as below: in which C 0 q denotes the initial covariance matrix of the proposal distribution q 1 and J is the Jacobian matrix, which is defined as below: The partial derivatives ∂y i (p) ∂p j (i = 1, 2, . .., m; j = 1, 2, . .., n) are calculated at the optimal solution of the Nelder-Mead simplex method.Because these derivatives cannot be solved using analytical methods, we approximated them using ∂y i (p) ∂p j ≈ (18) y i p 1 , . .., p j + p j , . .., p m − y i (p 1 , . .., p j , . .., p m ) p j , in which p j was set at 5 % of the parameter p j .
3 Results and discussion

Synthetic soil column description and boundary conditions
To test the value of the developed joint inversion approach under a range of conditions and assumptions, we performed several synthetic case studies using the numerical soil column illustrated in Fig. 3.The synthetic column was developed to mimic typical soil and petrophysical properties associated with a high-centered polygon at an intensive study transect (NGEE Arctic, Barrow, Alaska) (Fig. 4).The transect is 35 m in length and covers three typical topography types in Barrow, namely, high-centered (HCP), flat-centered (FCP) and low-centered polygon (LCP).The thawing occurs during the growing season, lasting from the beginning of June to the end of September.In the growing season We assumed that the vertical profiles of soil properties are constructed by interpolating their corresponding values at z = 0.1, 0.3, 0.6 and 1 m.For scenarios 1-8, the soil porosity was fixed ( = 0.9, 0.5, 0.5 and 0.8 for z = 0.1, 0.3, 0.6 and 1 m, respectively).For scenario 9, the soil porosity was calculated from the OC and sand content ( = 0.86, 0.67, 0.57 and 0.47 for z = 0.1, 0.3, 0,6 and 1 m, respectively).
while the LCP is fully saturated, the HCP is relatively dry and unsaturated.The bottom of the thaw layer at the end of the growing season is located at about 0.3 and 0.5 m depth at the center of HCP and LCP, respectively.ERT measurements were performed along the transect daily using the Wenner-Schlumberger configuration with an electrode spacing of 0.5 m.Other measurements and conditions useful for our synthetic studies -including soil temperature, soil moisture, thaw depth, snow dynamics and climate conditionswere also measured (Dafflon et al., 2017).These data have been used here to develop conceptual models and synthetic columns, while they will be used for real application of the joint inversion scheme in a subsequent study.Soil properties and petrophysical information used for the synthetic studies are provided in Table 1.The "true" soil properties are based on the core sample analysis at the Barrow, AK site (Baptiste Dafflon, personal communication, 2016), and the "true" petrophysical parameters were obtained from Minsley et al. (2015).It is worth noting that soil is rep- resented in the CLM as a mixture of OC, sand and clay.As such, in order to estimate the soil mixture, it was sufficient for us to consider OC and sand content (in sand-clay mineral mixture) only.
We assumed that the vertical profiles of soil properties (porosity, OC and sand content) were constructed by interpolating their corresponding values at four depths z k = 0.15, 0.3, 0.6 and 1 m as below: , where f z are the soil properties at depth z and f k are the soil properties at the corresponding depths z k = 0.15, 0.3, 0.6 and 1 m.These depths were chosen to represent the vertical variations of OC content and soil porosity in the core samples collected at the NGEE Arctic Barrow, Alaska, site.We synthetically explored nine scenarios using the newly developed inversion procedure (Table 2).The purposes of these scenarios are as follows: Scenarios 1 and 2: Evaluate the effect of measurement errors on uncertainties of soil OC estimates (using electrical resistivity data as an example).
Scenarios 2, 3, 4, 5, 6 and 7: Investigate the improvement in OC estimation gained by joint inversion of multiple hydrological, thermal and geophysical datasets compared with inversion of each single dataset.
Scenarios 7 and 8: Study how the parameter estimates and their associated uncertainties change if, in addition to OC content, sand content and petrophysical parameters are unknown.
Scenario 8 and 9: Explore the effect of soil porosity on the parameter estimation by comparing two cases: (1) soil porosity profile is fixed and independent from the soil OC and sand content, and (2) soil porosity is defined as a function of OC and sand content.
Scenario 8: Analyze the uncertainties, non-uniqueness, correlation and convergence of the inverse problem, and evaluate the impact of parameter uncertainty on prediction of hydrological-thermal dynamics.
For all scenarios, we used daily time step meteorological forcing data (including air temperature, wind speed, shortwave and long-wave radiation, and precipitation) collected at the Barrow site over a year period from 1 January to 31 December 2013, which includes a time period over which some of the soil and electrical datasets were also collected at the NGEE Arctic site.The plant functional type information was obtained from the CLM database for the Arctic region.The general approach that we followed to perform all synthetic scenarios is presented in Fig. 5.
In order to account for the measurement errors, we assumed that the error distribution was Gaussian and added error to synthetic data to obtain "noisy" synthetic data (hereafter referred to as observation data) (Table 2).We set the standard deviation of ERT measurement error to 2 % of synthetic data for scenario 1 (low measurement error) and to 5 % for the other scenarios.We used a standard deviation of measurement errors of 0.5 • C for soil temperature and 0.08 for soil liquid water content.The standard deviation of liquid water content was set to be relatively high because we observed that the associated error measurements for this variable at Barrow were quite high.The standard deviation of temperature measurement was set higher than what is generally expected for such measurements, while it also includes some error associated with relating measurements to precise depths.Observation data for inversion include (1) apparent resistivity data at the seven most important time points during the year, which include events such as thawing (day 163), summer growing season (days 185, 199 and 234), freeze-up (day 266) and frozen winter (days 292 and 312); (2) soil temperature data at z = 0.004, 0.16, 0.8, 1 and 2.4 m from day 49 to 365, which are the most varying; and (3) liquid water content at depths 0.004, 0.05, 0.11 and 0.2 m during the summer growing season (days 159 to 259).Liquid water in the winter season was not considered because it approximately equals zero and, therefore, does not contain any information for inversion.Meanwhile, soil temperature in the winter season still exhibits a large spatiotemporal variation, so we used the temperature data from both winter and summer for inversion.
For inversion, ranges were provided for unknown soil and petrophysical parameters based on Hubbard et al. (2013) and Dafflon et al. (2017) (Table 2).To minimize non-uniqueness in the inversion procedure, we ignored the small OC content at 1 m and the small sand content at 0.1 m.For scenarios 1-7, we estimated OC content at z = 0.1, 0.3 and 0.6 m (three parameters).For scenarios 8 and 9, we estimated OC content at z = 0.1, 0.3 and 0.6 m; sand content at z = 0.3, 0.6 and 1 m; and petrophysical parameters m and α (eight parameters).We assumed that there is no a priori information on the estimated parameters.As a result, the a priori distributions of OC and sand content were uniformly distributed within their parameter ranges.

Simulation results
In order to estimate the a posteriori PDF of OC and sand content as well as petrophysical parameters, we generate 8000 samples for scenarios 1-7 and 15 000 samples for scenarios 8 and 9.The number of samples in scenarios 8 and 9 is larger because there is a greater number of estimated parameters in these scenarios.We selected the last 5000 samples having a Geweke's score less than 0.4 to construct the PDFs of these parameters.Their best estimates and associated uncertainties are, respectively, represented by the means and standard deviations of the samples and summarized in www.the-cryosphere.net/11/2089/2017/The Cryosphere, 11, 2089-2109, 2017 Table 2. Discussion and comparison of the scenarios are presented below.

Effect of measurement error on parameter uncertainty
The influence of measurement error on the parameter uncertainties was considered by comparing scenarios 1 and 2 using apparent resistivity as an example.Scenario 1 assumed that the standard deviation of measurement error is 2 % of synthetic apparent resistivity data (small measurement error), while this value for scenario 2 is 5 % (large measurement error).For these two scenarios, we estimated the OC content at z = 0.1, 0.3 and 0.6 m. Figure 6 shows the probability functions of the OC at these depths.The figures indicate that the uncertainties of the estimated OC content at z = 0.1 and 0.3 m are considerably higher when the measurement error is larger.As shown in Table 2, when the measurement error of apparent resistivity increases from 2 to 5 %, the standard deviation of the a posteriori OC samples increases by a factor of 3 from 0.2 to 0.6 at z = 0.1 m and by more than a factor of 2 from 2.6 to 5.5 at z = 0.2 m.At z = 0.6 m, the OC content cannot be reliably obtained by both scenarios.
In order to investigate the non-uniqueness problem and the correlation between parameters, we estimated the misfit (sum of square of absolute differences) between the synthetic and sampled apparent resistivity data as a function of the OC content at z = 0.1, 0.3 and 0.6 m for scenario 2 (Fig. 7).While the OC at z = 0.1 m is well identified, the misfit negligibly changes when the OC content at z = 0.6 m varies from 20 to 50 %.This indicates that the apparent resistivity data are insensitive to OC content at z = 0.6 m.This is reasonable, because this depth is within the permafrost (see Fig. 13), where temperature insignificantly changes over time.Results could be different if the permafrost were deeper.Indeed, OC content at depths inside the active layer where a freeze-thaw process occurs is expected to be better resolved because of the stronger temporal changes in properties used in the inversion.Figure 7 also shows that there is a negative correlation between the OC content at z = 0.3 and 0.6 m, which increases the uncertainties of OC estimates at both depths.

Influence of joint inversion of multiple data on parameter uncertainty
The effectiveness of the joint inversion of multiple datasets on the OC content estimation (at z = 0.1, 0.3 and 0.6 m) was investigated by comparing results obtained from six scenarios that used (1) single apparent resistivity (scenario 2); (2) single temperature (scenario 3); (3) single liquid water content (scenario 4); (4) temperature and apparent resistivity data (scenario 5); (5) liquid water content and apparent resistivity (scenario 6); and (6) liquid water content, temperature and apparent resistivity data (scenario 7).Joint inversion of apparent resistivity with temperature and/or liquid water con- tent data significantly reduces the uncertainties of OC content at z = 0.1 and 0.3 m (Fig. 8).For example, compared to using a single temperature dataset, the uncertainty of OC content (the standard deviation of the final Markov chain of OC content) reduces from 0.4 to 0.2 at z = 0.1 m and from 4.6 to 1.9 at z = 0.3 m when jointly using temperature, liquid water content and apparent resistivity datasets.Finally, we found that, even when all "observation" data are used, there is no improvement in the OC content estimate at z = 0.6 m.These synthetic experiments suggest that, given this depth is located within the permafrost (see Fig. 13), the apparent resistivity, liquid water content and temperature data are insensitive to OC content.This is because within the permafrost the soil temperature and ice-liquid water content exhibit much smaller variations than in active layer, in both time and space.

Effect of mineral content and petrophysical parameters
In scenario 8, in addition to the OC content, we assumed that the sand content and petrophysical parameters m and α are unknown and estimated these parameters using the apparent resistivity, temperature and liquid water content data.Similar to the previous scenarios, the OC content at z = 0.1 and 0.3 m was obtained with small uncertainties (σ OC (z = 0.1 m) = 0.3, σ OC (z = 0.3 m) = 2) (Fig. 9).The sand content at z = 0.3 and 1 m was also well estimated with uncertainties of 2.4 and 2, respectively.It is worth noting that, regardless of deep location, the sand content at 1 m is relatively well determined because at this depth the sand-clay mineral (92 %) dominates the OC content (8 %), and therefore the hydrological-thermal data are relatively sensitive to this parameter.By contrast, the OC and sand content at z =0.6 m is unidentifiable with uncertainties up to 6.5 and 8.2, respectively.Finally, both of the petrophysical parameters m and α are well estimated, while the parameter α has lower uncertainty.This implies that α is more sensitive to the apparent resistivity than to m.The pairwise relationships between estimated parameters (Fig. 10) indicate that the OC content at 0.1 m and petrophysical parameter α are the most reliably estimated parameters, followed by the OC content at z = 0.3 m, sand content at z = 0.3 and 1 m, and cementation index m.As for the correlation between parameters, the figure reveals that there is a strong positive correlation between the sand and OC content at z = 0.6 m with a correlation coefficient of 0.86.This correlation and the insensitivity of the observations with their variations are two main reasons for the non-uniqueness of these two parameters.The pairs of m-α and the OC z = 0.1 msand z = 0.3 m are also highly correlated, with correlation coefficients of 0.84 and 0.70, respectively.

Effect of porosity dependence on OC and mineral content
In this section, we evaluate how the parameter uncertainties change when the porosity is determined as a function of the OC and mineral content by comparing scenarios 8 and 9.
While the soil porosity in scenario 8 was fixed and independent from the OC and sand content, it was calculated from the OC and sand content in scenario 9 as shown in Eqs. ( 2) and (3).Compared to scenario 8, all uncertainties of sand and OC content in scenario 9 are smaller (Fig. 11).In particular, the uncertainties of these parameters at z = 0.6 m significantly decrease from 6.5 to 3.8 (for OC content) and from 8.2 to 1.8 (for sand content).This can be explained by the fact that, in addition to thermal parameters (thermal conductivity and heat capacity), the OC and sand content in scenario 9 controls the soil porosity, which also influences the subsurface hydrological-thermal dynamics (see Fig. 2).As a result, the temperature, liquid water and apparent resistivity data in this scenario are more sensitive to variations of OC and sand content than those in scenario 8. Consequently, these parameters are more identifiable.By contrast, the uncertainty of petrophysical parameters α and m considerably increases from 0.002 to 0.022 (for α) and from 0.042 to 0.066 (for m).This is because, while it was fixed in scenario 8, the soil porosity depends on the OC and sand content in scenario 9. Therefore, the soil porosity in scenario 9 is also uncertain due to the uncertainties of the OC and sand content.Because the soil porosity, α and m are closely correlated (see Eq. 4), the uncertainty of soil porosity causes higher uncertainties of α and m.

Uncertainty propagation from parameters to the hydrological-thermal and thaw layer thickness prediction
In this section, we evaluate the impact of parameter uncertainties on the prediction of hydrological-thermal dynamics.
A posteriori samples of the OC, sand content, and petrophysical parameters m and α of scenario 8 were used for this analysis.The synthetic and estimated soil temperature at z = 0.004, 0.16 and 0.99 m and the liquid water content at z = 0.004, 0.05 and 0.11 m are compared (Fig. 12).The uncertainties of these predictions are represented by grey regions with a confidence interval of 95 %.The figure indicates that the synthetic and estimated soil temperature and liquid water content agree well with each other.However, the uncertainty of the soil temperature prediction is much smaller than that of the liquid water content prediction.The average confidence intervals over the simulation period of the soil temperature prediction at z = 0.004, 0.16 and 0.99 m are 2.3, 2.3 and 1.7 % of the "observation", respectively, while these values for the soil liquid water content prediction at z = 0.004, 0.05 and 0.11 m are 28.7,16.1 and 12.9 %, respectively.These dif-    ferences can be explained by the high sensitivity to the OC and sand content and the larger measurement errors of liquid water content than of soil temperature.The synthetic and estimated thaw depth using results obtained from scenario 8 (Fig. 13) show that soil water thaws around the middle of June and freezes again around the middle of September.The thaw depth varies from 0.2 to 0.42 m.These results are compatible with our field survey data in Barrow (Dafflon et al., 2017), indicating that, although this is a synthetic study, its simulation is relatively compatible with the Arctic tundra field measurements.As for the influence of parameter uncertainties on the thaw depth estimation, we observed that the parameter uncertainties only cause thaw depth variations during the warmest period of the year (beginning of August to middle of September).During other times of the year, the thaw depths corresponding to different sets of parameters are quite similar.
The comparison between synthetic and predicted apparent resistivity data (Fig. 14) shows that there is a very good agreement between them with no bias, which implies that our inversion scheme converges to the lowest misfit region.The confidence ranges corresponding to a level of 95 % vary from 1.4 to 9.4 % of the "observation" resistivity, which is suitable with the relative measurement error of 5 %.

Summary and conclusions
In this study, we developed and tested a surface-subsurface coupled hydrogeophysical inversion approach to estimate OC content and its influence on hydrological-thermal behavior under Arctic freeze-thaw conditions.In our inversion scheme, the CLM model serves as a forward model to simulate the land surface energy balance and surface-subsurface hydrological-thermal processes.The new scheme can jointly use different types of data for the inversion, including electrical resistivity data.The dependence of soil electrical resistivity on temperature and ice-liquid water content is explicitly accounted for within the inversion.
We developed an advanced optimization technique that combines the deterministic and stochastic optimization algorithms to obtain soil and petrophysical parameters and their associated uncertainties.The stochastic optimization estimated the a posteriori distribution of model parameters by using the Bayesian inference and adaptive MCMC algo-rithm DRAM.Meanwhile, the deterministic optimization algorithm was used to approximate the starting set of model parameters and the initial covariance matrix of the proposal distribution for the stochastic optimization, which helps to more quickly converge to the parameter a posteriori distribution.
We tested the inversion scheme using multiple synthetic experiments in a 1-D soil column representative of the Arctic tundra, where surface-subsurface hydrological and thermal regimes co-interact and are influenced by soil OC and mineral content.The obtained results show that the new inversion approach reproduced the synthetic data well in all experiments.The shallow (upper 0.3 m) active-layer OC and sand content and the petrophysical parameters can be reliably obtained using soil temperature, soil liquid/ice water content and ERT data.When the soil porosity is fixed, the uncertainties of OC and sand content are very high in the permafrost section (0.6 m), even when soil temperature, liquid water saturation and apparent resistivity data are jointly used in the estimation procedure.This suggests that, when the porosity is fixed, the inversion approach is unable to significantly improve the estimation of OC within the permafrost, due to the small magnitude of temporal variation of both temperature and soil moisture in that section.However, if the soil porosity is considered as a function of OC and sand content, the permafrost parameters can be reliably obtained because the variation of porosity with OC and sand content increases the sensitivity to ice-liquid water and temperature.Examining the relationship between measurement errors and parameter uncertainties, we found that the uncertainties of estimated parameters increase with increasing measurement error.We also explored the improvement in parameter estimation when jointly using multiple data for the inversion.Compared to single dataset inversion (temperature, soil moisture or electrical resistivity), joint inversion significantly reduces the uncertainties of estimated parameters, especially at 0.3 m depth.Finally, we quantified the influence of parameter uncertainties on the prediction of hydrologicalthermal and thaw depth dynamics.The obtained results show that the soil liquid water content prediction is more uncertain than the soil temperature and apparent resistivity predictions, due to its large measurement error.The uncertainties in OC and sand content have an impact on the thaw depth estimation only during the warmest months of the year (August and September).
This study developed and tested a novel approach to estimating soil OC content using inverse modeling that can incorporate diverse hydrological, thermal and ERT datasets.In addition, the study also permitted exploration of surfacesubsurface hydrological-thermal dynamics and spatiotemporal variations associated with freeze-thaw transitions.Given the importance of characterizing OC as part of ecosystem and climate studies, the typical challenges associated with collecting and analyzing "sufficient" core data to characterize the vertical and horizontal variability of OC associated   with a field study site, and the increasing use of electrical resistivity data to characterize vertical, horizontal and temporal variability in shallow systems, the new inversion approach offers significant potential for improved characterization of OC over field-relevant conditions and scales.It also offers significant potential for improving our understanding of hydrological-thermal behavior of naturally heterogeneous permafrost systems.The successful validation of this approach using 1-D synthetic studies provides a foundation for extending it to 2-D and applying it to real field data, which is currently underway.
In this study, we concentrated on the indirect impact of the OC content on water electrical resistivity via soil water and temperature.Recent studies indicated that the soil OC content largely influences ionic mobility and, therefore, changes the polarization and relaxation time of soil response to the applied current, which can be measured by spectral induced polarization (SIP) (e.g., Schwartz and Furman, 2015).As a result, our future study will explore the possibility of integrating SIP measurements into our coupled hydrological-thermal-geophysical inversion scheme.In that case, the OC content is linked to SIP measurements both by its hydrological-thermal and by its electrical polarization properties.Hauck et al. (2011) indicated that combination of ERT and seismic measurements can improve the estimation of ice and liquid water.We will integrate this approach into coupled hydrogeophysical inversion to better constrain the inversion and reduce the non-uniqueness of parameter estimation.
With advancements in data acquisition, the surfacesubsurface hydrological-thermal dynamics now can be monitored in real time at high temporal resolution using multiple above-and below-ground measurements including geophysical techniques.Our next step is to expand the inversion scheme so that it can assimilate these data into hydrologicalthermal models to improve the model prediction in real time.

Figure 1 .
Figure 1.(a) Forward coupled hydrological-thermal-geophysical model that considers soil liquid/ice water content, temperature and apparent resistivity.(b) The two-stage inversion scheme combines deterministic and stochastic optimization algorithms to estimate the PDFs of desired model parameters (p), which include the OC content (scenarios 1-9), sand content (scenarios 8 and 9) and petrophysical parameters (scenarios 8 and 9) (see Table2).The scheme permits flexibly using single or multiple types of data for inversion.The forward coupled hydrological-thermal-geophysical model (a) is iteratively executed in both deterministic and stochastic inversion stages.

Figure 2 .
Figure 2. Soil thermal conductivity (a, d), heat capacity (b, e) and thermal diffusivity (c, f) as a function of the liquid water saturation, OC and sand content.The calculation for this figure was based on equations presented in Appendix A. The soil porosity was fixed at 0.7 for (a)-(c) and determined as a function of OC and sand content (Eqs. 2 and 3) for (d)-(f).(a) Thermal conductivity; porosity φ = 0.7.(b) Volumetric heat capacity; porosity φ = 0.7.(c) Thermal diffusivity; porosity φ = 0.7.(d) Thermal conductivity; porosity is a function of OC and sand content.(e) Volumetric heat capacity; porosity is a function of OC and sand content.(f) Thermal diffusivity; porosity is a function of OC and sand content.

Figure 3 .
Figure3.Twenty-seven synthetic soil layers and soil properties (OC, sand content and porosity) for simulating hydrological and thermal dynamics.The five bottom bedrock layers are not shown in this figure.We assumed that the vertical profiles of soil properties are constructed by interpolating their corresponding values at z = 0.1, 0.3, 0.6 and 1 m.For scenarios 1-8, the soil porosity was fixed ( = 0.9, 0.5, 0.5 and 0.8 for z = 0.1, 0.3, 0.6 and 1 m, respectively).For scenario 9, the soil porosity was calculated from the OC and sand content ( = 0.86, 0.67, 0.57 and 0.47 for z = 0.1, 0.3, 0,6 and 1 m, respectively).

Figure 4 .
Figure 4. (a) Image of the intensive ERT transect (dash line) and pole-mounted cameras, which monitor the land surface variability of the whole transect.(b) Aerial view of the ERT transect (dashed line), which covers different types of polygons, namely, high-centered polygon (HCP), flat-centered polygon (FCP) and low-centered polygon (LCP).(c) An example of inversion of ERT data measured in August 2013.The top layer with low resistivity represents the active layer.The middle layer with high resistivity and the underlying less resistive layer correspond to permafrost and saline permafrost, respectively (see Dafflon et al., 2017, for more details).

Figure 5 .
Figure 5.General procedure used to perform the synthetic case studies.

Figure 6 .
Figure 6.The a posteriori probability of the soil OC content at z = 0.1, 0.3 and 0.6 m obtained by inverting apparent resistivity data with relative measurement error of 2 % (a) and 5 % (b).The red lines represent the "true" OC content.

Figure 7 .
Figure 7. Misfit (sum of square of absolute difference) between synthetic observations and MCMC sampling apparent resistivity data as a function of soil OC content at z = 0.1, 0.3 and 0.6 m for scenario 2.

Figure 9 .
Figure 9.The a posteriori probability of soil OC content at z = 0.1, 0.3 and 0.6 m; sand content at z = 0.3, 0.6 and 1 m; and petrophysical parameters m and α for scenario 8.The sand content is the fraction of sand in the sand-clay mineral mixture.Soil porosity is fixed and independent from the OC and sand content.The red line represents the "true" parameter values.

Figure 10 .
Figure 10.Pairwise relationships between estimated parameters.The calculation was based on 3000 MCMC samples of scenario 8.

Figure 11 .
Figure 11.The a posteriori probability of soil OC content at z = 0.1, 0.3 and 0.6 m; sand content at z = 0.3, 0.6 and 1 m; and petrophysical parameters m and α for scenario 9.The sand content is the fraction of sand in the sand-clay mineral mixture.Soil porosity is determined as a function of soil OC and sand content in the CLM.The red line represents the "true" parameter values.

Figure 12 .
Figure 12.Comparison of "observation" and predicted soil temperature at z = 0.004, 0.16 and 0.99 m (a) and liquid water content at z = 0.004, 0.05 and 0.11 m (b).The blue line denotes the synthetic data.The grey region represents the 95 % confidence interval calculated from the a posteriori MCMC samples of scenario 8.The red line represents the mean of samples.

Figure 13 .
Figure 13.Comparison between estimated and synthetic thaw depth over a year for scenario 8.The blue and red lines, respectively, represent the synthetic and estimated thaw depth.The grey region shows the confidence interval with a level of 95 %.

Figure 14 .
Figure 14.Comparison between "observation" and predicted apparent resistivity.The red line denotes the 1 : 1 line.The vertical error bar on the blue symbols represents the confidence interval of the predicted apparent resistivity with a confidence level of 95 %.The comparison was based on a posteriori samples of scenario 8.