2.1 ATS–BERT model

To set up the synthetic model, we used digital elevation data of a transect through ice-wedge polygonal tundra at the Barrow Environmental Observatory (BEO) at Utqiaġvik, Alaska (Fig. 2). Our study includes an 11 m section covering a single polygon with an ice wedge on each side. In this study we do not explicitly assign ice properties for the ice wedges. Instead, we model bulk porosities and effective thermal conductivities that can be associated with peat and mineral layers of the entire transect.

In Fig. 2a, we present the computational mesh representing the cross section of the polygonal tundra that ATS is run on. The thickness of the peat layer corresponds to observations at the site, with a thick peat layer on the sides (troughs) and a thinner layer in the middle of the low-centered polygon. A mineral layer was assigned below the peat layer across the transect. We initially designated six synthetic direct temperature and soil moisture measurement locations within the active layer area, the maximum thaw layer from the ground surface to the top of the permafrost, similar to the sensor setup at the site (Dafflon et al., 2017). The average active layer depth is about 38 cm, as it can be seen from the ground temperatures simulated for the synthetic model run with actual meteorological data in Fig. 2b. The linear white region in Fig. 2b indicates the bottom of the active layer within the transect (0 ^∘C). Then we added four more synthetic direct measurement locations below the active layer to evaluate the effect of their inclusion on PE accuracy and robustness. All observation locations are represented as stars in Fig. 2a corresponding to the locations of the collected daily averaged temperature and soil moisture time series. The temperature and soil moisture time series were recorded at depths of 5, 20, 60, and 80 cm below the surface.

Figure 2The (vertically exaggerated) 2-D transect used by the ATS model. (a) The unstructured mesh where green represents the peat layer and brown represents the mineral soil layer. Black stars represent the six sensors recording temperature and soil moisture content within the active layer. Red stars represent the four sensors recording temperature and soil moisture content below the active layer. (b) Ground temperature distribution simulated by the ATS model, corresponding to the time of maximum active layer depth. Here the depth of the active layer corresponds to the distance above the white linear feature (i.e., 0 ^∘C) dividing the thawed and frozen regions of the ground. The light blue dots represent the location of the electrodes in this setup.

Download

The setup of the ATS model followed a standard procedure described in several studies (Atchley et al., 2015; Painter et al., 2016; Jafarov et al., 2018). Typically, we set up the model in several steps: (1) initialization of the water table, (2) introduction of the energy equation to establish antecedent permafrost, and (3) spin-up of the model with simplified and actual meteorological data from the BEO station. We spun up the model until the active layer achieved cyclical equilibrium. The overall depth of the modeling domain is 50 m. We set the bottom boundary to a constant temperature of T=263.55 K and set zero heat and zero mass flux boundary conditions on the vertical sides. A seepage face was imposed at 4 cm below the surface on each side of the domain to allow drainage to the trough network to prevent water from pooling at the surface, as is typical of partially degraded polygonal ground (Liljedahl et al., 2016). We use two types of meteorological datasets as surface boundary condition drivers for the ATS model: simplified (sinusoidal air temperature, constant precipitation, and constant radiative forcing) and actual weather data from the BEO site. The actual meteorological data were collected starting on 1 January 2015 and include air temperature, rain and snow precipitation, humidity, long- and shortwave radiation, and wind speed.

We created a synthetic “truth” by designating porosities and soil grain thermal conductivities {ϕ,k} of peat and mineral soil as parameters in the forward model. The resulting temperature (T), liquid and ice saturations (s_l, s_i), and apparent resistivities (ρ_a) were collected as the true state. Critical for these simulations is the calculation of the thermal conductivities of the bulk soil; calculated in ATS using Kersten numbers to interpolate between saturated frozen, saturated unfrozen, and fully dry states (Painter et al., 2016) where the thermal conductivities of each end-member state are determined by the thermal conductivity of the components (soil grains, air, water, or liquid) weighted by the relative abundance of each component in the cell (Johansen, 1977; Peters-Lidard et al., 1998; Atchley et al., 2015). Thermal conductivities of water, ice, and air are considered constant, leaving soil grain thermal conductivity as the remaining parameter to be estimated. The equation to calculate saturated, frozen thermal conductivity (κ_sat,f) has the following form:

where κ_sat,uf, κ_i, κ_w are thermal conductivities for saturated unfrozen, ice, and liquid water, respectively; and ϕ is porosity.

The freezing characteristic curve is thermodynamically derived using a Clapeyron relation and the unfrozen water retention curve, as described in Painter and Karra (2014) and Painter et al. (2016). In Fig. 3 we present liquid and ice saturations for one realization of the model for winter (January) and summer (August) times of the year. The ice saturation is high below the active layer all year long and lowest within the active layer in the summer. The peat layer holds more water; therefore, ice concentration is higher than in the mineral layer in the winter. The liquid saturation plot shows that, by the end of the summer, the peat layer is drier than the mineral layer.

We sequentially couple the ATS and BERT numerical models via petrophysical relationships used by Tran et al. (2017) and based on Archie (1942) and Minsley et al. (2015). In that approach, the electrical resistivity (ρ) is determined as a function of soil characteristics, temperature, porosity, liquid water saturation, fluid conductivity, and ice content:

where σ_w is the fluid electrical conductivity, σ_s is soil/sediment electrical conduction, n is a saturation index, d is a cementation index, and c is a temperature compensation factor accounting for deviations from T=25 ^∘C.

The ice content is linked to water content through the liquid water saturation and to σ_w, which is influenced by the concentration of Na⁺ and Cl⁻ ions and the ice–liquid fraction. Here σ_w has the following form:

where F_c is Faraday's constant, β_i and z_i are the ionic mobility and valence respectively, C_i is the concentration of the ith ion, and α is the factor influencing how the liquid water salinity increases when the fractions of liquid in ice–liquid water $S_{f_{\text{w}}}$ decrease. $S_{f_{\text{w}}}$ is defined as

Both s_l and s_i are simulated by ATS. Note that ϕ in Eq. (2) is an estimated parameter (see Fig. 1). In this study we assume that the n, d, σ_s, α, F_c, ${\mathit{\beta }}_{{\mathrm{Na}}^{+}}$, ${\mathit{\beta }}_{{\mathrm{Cl}}^{-}}$, $C_{{\mathrm{Na}}^{+}}$, and $C_{{\mathrm{Cl}}^{-}}$ parameters used in Eqs. (2) and (3) are known (see Tran et al., 2017) and focus on the robustness of the PE algorithm in estimating porosity and thermal conductivity.

Figure 3The 2-D transects of daily snapshots simulated by the ATS model. The rows from top to bottom correspond to ice saturation (a, b) and liquid saturation (c, d). The columns from left to right correspond to a day in the middle of the winter and the day of maximum active layer depth in summer.

Download

The 2-D resistivity data inferred from ATS simulations and petrophysical relationships get passed to BERT, which simulates resistances that are then converted to the apparent resistivities (ρ_a). The ρ_a values correspond to an acquisition along an 11 m long transect using a 0.5 m electrode spacing and a Schlumberger configuration with a total of 138 measurements (see Fig. 2b). This configuration implies that the measurements are mostly sensitive to the electrical resistivity in the top few meters.

Since BERT and ATS operate on different unstructured meshes, we wrote a function that interpolates the values between the two meshes. Note that the ATS mesh is 50 m deep. We calculate ρ by using corresponding outputs from the ATS model and the petrophysical relationships, and we then interpolated these values on a mesh defined in BERT and adapted to the acquisition geometry. BERT's mesh consists of a finely resolved mesh (11 m wide by 4.5 m deep) embedded within a coarser outer mesh that is about 120 m wide and 85 m deep. We link hydrological variables with electrical resistivities in the fine mesh. The coarse mesh is used to reduce the effect of boundaries. It extends until the change in the electrical resistivity between two neighboring cells is negligible.

2.2 Parameter estimation using PEST

To test if the known soil properties can be recovered by the PE approach, we start with randomly selected initial parameter guesses. We use a Latin hypercube sampling method to generate random initial guesses of porosity and thermal conductivities around the synthetic truth (McKay et al., 1979). Each parameter combination includes four parameters: porosity and thermal conductivity for peat and mineral soil layers. These parameters were chosen due to their strong controls on both hydrologic and thermal states (Atchley et al., 2015; Nicolsky et al., 2009). The rest of the hydrothermal properties are kept fixed.

The inverse approach involves the minimization of a cost function expressed as the sum of squared differences between simulated values and synthetic measurements using the Levenberg–Marquardt (LM) algorithm (Levenberg, 1944; Marquardt 1963) implemented in the PEST software package (Doherty, 2001), which was used to handle all parameter estimation runs.

To estimate soil porosities and thermal conductivities, we minimize the cost function (J), which includes calculated and synthetic T, s_l, and ρ_a in the following form:

where subscripts c and s correspond to calculated and synthetic states of the system; and w_T, w_s, and $w_{{\mathit{\rho }}_{\text{a}}}$ are the corresponding weights for the temperature, saturation, and apparent resistivity residuals. n_sens is the number of sensors, n_days is the number of days over which we collected the data, n_snap is the number of ρ_a snapshots, and n_meas is the number of ρ_a measurements during one snapshot. T_c and $s_{{\text{l}}_{\text{c}}}$ are time series from multiple sensors collected daily from the beginning of June until the end of September. ρ_a represents the apparent resistivity data snapshots taken at a certain day. The number of apparent resistivity snapshots depends on the particular case, varying from one to eight snapshots per year. The one-snapshot case corresponds to only one snapshot in the month of August, while the eight-snapshot case corresponds to a snapshot taken once per month from January until September. In addition, we tested the case where we collected eight daily ρ_a snapshots. This was done to compare how different time spacing would affect the estimated properties.

The weights were chosen in order to scale the contribution of each type of residual so that contributions to the cost function are evenly distributed across temperature, saturation, and apparent resistivity residuals. For example, saturation residuals are on the order of a few tenths, while apparent resistivity residuals can be tens of ohmmeters. The weights were selected based on evaluating the individual contributions to the cost function for each measurement type on an ensemble of simulations spanning the parameter ranges. The apparent resistivity residual weight ($w_{{\mathit{\rho }}_{\text{a}}}$) was set to one. The temperature and saturation residual weights (w_T and w_s) were then modified so that each measurement type component in the cost function had roughly equivalent magnitude over most of the parameter space. This resulted in weights of $w_{{\mathit{\rho }}_{\text{a}}}=\mathrm{1}$, $w_{\text{T}}=\sqrt{\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{3}}}$, and $w_{\text{s}}=\sqrt{\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{5}}}$.

If the cost function satisfies a minimum criterion or the maximum allowed number of iterations, which we chose to be equal to 25, is reached, the PE terminates. The porosities and thermal conductivities corresponding to the minimum of the cost function, i.e., the parameters associated with the best fit between simulated and synthetic values, are considered the estimated parameter values as

Here {ϕ,k} are the estimated porosities and thermal conductivities for peat and mineral soil.

Based on sensitivity analyses using simplified meteorological data, the cost function response surface was smooth and convex over the parameter ranges of interest. Therefore, we chose the LM approach because of its robust gradient-based optimization scheme that takes advantage of smooth convex response surfaces to quickly converge to minima.

Table 1Allowed range for the estimated parameters.

Download Print Version | Download XLSX

Figure 4Estimated properties from 30 inversions of the (S)3T3s_l1ρ_a case, where the “true” values are shown as a cross section of two dashed lines for the bulk porosities and effective thermal conductivities for peat and mineral soil layers. Yellow lines correspond to the paths taken by the LM algorithm. The white dots correspond to the estimated values. (a) Projection of the cost function with respect to porosities of the peat and mineral layer. (b) Projection of the cost function with respect to thermal conductivities of the peat and mineral layer. The color bar represents the cost function normalized by its maximum logarithmic value.

Download

2.3 Experiments

To build an understanding of the inverse framework, we start with a simple setup and then gradually add more complexity. First, we use simplified meteorological data where we assume that air temperatures change according to a sinusoidal function and all other terms are constants. Initially we start with three temperature and moisture content measurement locations within the peat layer (refer to Fig. 2a) and one ERT data snapshot. Then we increase the number of ERT data snapshots up to eight by adding snapshots once per month from January until August. Each ERT data snapshot calculated by BERT uses the set of daily averaged T and s_l simulated by ATS and petrophysical relations (Eqs. 2 and 3) which are varying over time. Then we increase the number of sensors to six and add noise to the simulated data. Introduction of noise allows us to evaluate the effect of measurement uncertainties that will be present in the actual application of the PE method. We added different levels of Gaussian noise to the synthetic measurements of T, s_l, and ρ_a in the following way: 1 % to T, 5 % to s_l, and 10 % to ρ_a. These levels of noise for the different types of measurements are based on published literature and our own experience (Wang et al., 2018; Dafflon et al., 2017). After that we substitute simplified meteorological data with actual data from the BEO site to evaluate our PE method under realistic ground surface boundary conditions. In this case we evaluate how much and what kind of data we need to robustly recover subsurface porosities and thermal conductivities. To do this we test the inclusion of individual types of measurements in the cost function (Eq. 3) as well as all possible combinations of measurement types. We used different soil property ranges for the simplified and actual meteorological data PE runs which are summarized in Table 1. This was done to ensure that PE is able to recover different sets of parameters and to test the consistency and effectiveness of the PE method. Finally, we compared the difference between estimated parameters for eight ERT data snapshots collected once a month versus once a day for 8 d. The notation and a description of each run for simplified and actual meteorological data are summarized in Table 2.

Table 2Description of all PE cases used in this study. The numbers before T and s_l correspond to the number of sensors used. The number before ρ_a corresponds to the number of apparent resistivity snapshots used. n stands for noise added to the synthetic measurements. (S) corresponds to runs driven by simplified meteorological data. (s) represents daily ρ_a snapshots.

Download Print Version | Download XLSX

3.1 Simplified meteorological data

To evaluate the PE method performance driven by simplified meteorological data, we ran PE experiments using 30 different random combinations of porosity and thermal conductivity values as the initial starting point. We used 30 PE samples of {ϕ,k} starting points in the first experiment ((S)3T3s_l1ρ_a) to illustrate the overall performance of the parameter estimation using a large number of samples. After that, we did only five PE runs for the simplified meteorological data and 10 for all other runs with actual meteorological data. For all figures after Fig. 4, for consistency and clarity, we show results for only five PE runs per case. It is important to note that the number of samples that one needs to run to ensure the robust convergence of the estimated parameters depends on the specifics of the corresponding case (i.e., experiment-specific). If most of the LM runs converge to the same set of parameters and have low cost function values, then that set of runs most likely corresponds to the actual {ϕ,k}. In Fig. 4, the red triangles represent initial guesses (parameter combinations), and the synthetic truth is indicated by the intersection of the two dotted lines. Yellow lines connecting red triangles with white crosses represent the path that the LM algorithm has taken from the initial guess to the estimated parameter combination (white crosses, Fig. 4). The yellow dots along the yellow lines indicate the location at each LM iteration.

Figure 4 indicates that the method is able to recover porosities more robustly than thermal conductivities, i.e., estimated porosities are similar to their true state. According to the liquid saturation plot in Fig. 3, liquid saturation of the mineral layer is quite dynamic and more saturated in comparison to the peat layer. Nevertheless, thermal conductivity of the mineral layer (k_m) corresponds to the highest mismatch. Three out of 30 inversions corresponding to k_m end up close to 1.4 W m⁻¹ K⁻¹ (the true value is 2 W m⁻¹ K⁻¹), suggesting those values do not correspond to the truth, since most of the estimated values (27 cases) are concentrated around the intersection of the dotted lines. The response surface for the corresponding cost function (Eq. 5) lies hereby in a flat, low-gradient region. The projections of the cost function response surfaces corresponding to porosities (Fig. 4a) have a better defined minimum, as opposed to projections of the cost function response surfaces corresponding to thermal conductivities (Fig. 4b), indicating the nonuniqueness of the estimated parameters. For this experiment, we used time series of T and s_l only from the first three near-surface sensors (Fig. 2a). All of these three sensors are located in the peat layer, suggesting that using just near-surface sensors only from one upper layer might not be enough to recover the deeper layer thermal conductivity.

To illustrate the effect of noise on the robustness of the estimated parameters we used cases with six near-surface sensors (6T and 6s_l) and a varying number of ERT snapshots driven with simplified meteorological data. Similarly to the (S)3T3s_l1ρ_a case, (S)6T6s_l1ρ_a shows good convergence for porosities and poor convergence for thermal conductivities with an averaged error of 0.1 W m⁻¹ K⁻¹ (Fig. 5a, b). Adding noise to the (S)6T6s_l1ρ_a+n case slightly worsens the estimated porosity values and significantly worsens k_m with root-mean-squared error (RMSE) rising from 10 % to more than 50 % (Fig. 5c, d). Figure 5e and f show that increasing the number of ERT snapshots from one to eight per year (i.e., collected once per month from January until September) improves k_m estimates, allowing better convergence for four out of five samples to the synthetic truth. If we compare all three cases in Fig. 5 on how well they are able to estimate k_m, it is clear from Fig. 5d that for the case (S)6T6s_l1ρ_a+n none of the k_m's were correctly estimated, whereas significantly improved k_m values were found by increasing the number of monthly ERT snapshots (Fig. 5f). Moreover, all except one estimated value showed a better match with its true value than the (S)6T6s_l1ρ_a case without any added noise.

Figure 5Estimated properties from five inversions of the three different cases, where the true values are shown as a cross section of the two dashed lines for the bulk porosities and effective thermal conductivities for peat and mineral soil layers. Yellow lines correspond to the paths taken by the LM algorithm. The white dots correspond to the estimated values. The rows from top to bottom correspond to cases (S)6T6s_l1ρ_a (a, b), (S)6T6s_l1ρ_a+n (c, d), and (S)6T6s_l8ρ_a+n (e, f) respectively. The columns from left to right correspond to the projection of the cost function with respect to porosities and thermal conductivities. The color bar represents the cost function normalized by its maximum logarithmic value.

Download

Figure 6Five matrix tables presenting fitness metrics between synthetic model values and values obtained by the parameter estimation method using simplified meteorological data. Matrix tables from left to right correspond to the normalized root-mean-squared errors for temperatures, liquid water saturations, and apparent resistivities as well as to the normalized Euclidian distances between synthetic (true) and estimated values of porosities and thermal conductivities. Each matrix value was normalized by dividing it by the matrix maximum value. The normalized values are indicated by tildes.

Download

In Fig. 6, we summarize results of the five PE runs for each of the first six cases corresponding to simplified meteorological data listed in Table 2. The first three matrix tables correspond to the normalized RMSE values for each measurement type (ΔT, Δs_l, and Δρ_a). The last two matrix tables correspond to the normalized Euclidian distances between the synthetic truth and estimated parameter values of δϕ and δk. We normalized the values in each matrix by dividing by the maximum value from the corresponding matrix. The normalized values are marked with tildes and range from 0 to 1, where values closer to 0 correspond to a better match and values closer to 1 correspond to a worse match. As shown above, the method is able to accurately estimate both peat and mineral soil porosities as well as peat layer thermal conductivity (k_p), but it cannot always accurately estimate k_m. There is not much difference between cases (S)3T3s_l1ρ_a and (S)6T6s_l1ρ_a except for a slight improvement in k_m, suggesting that the small vertical distance (10 cm) between sensors 1 and 4, 2 and 5, and 3 and 6 could limit the recorded data variability, leading to difficulties in the estimation of the k_m parameter. Since all six sensors lie within the active layer, we added additional sensors below the active layer in the later experiments (red stars in Fig. 2a). The ϕ and k matrix tables show that increasing the number of monthly ERT snapshots consistently improves the estimates of ϕ and k. This suggests that increasing the number of monthly ERT snapshots can lead to improved convexity of the cost function (Eq. 5).

3.2 Meteorological data from the Utqiaġvik (Barrow) site 2015

After testing the PE method for the simplified meteorological data, we applied measured meteorological data from the BEO site for the year 2015. To better understand the importance of each measurement type and their combinations within the developed PE algorithm, we tested all of the scenarios corresponding to the “actual meteorological data” column from Table 2. The results of these runs are summarized in the colored matrix tables in Fig. 7. Since there are more than twice the number of actual meteorological cases than simplified meteorological cases, it is difficult to analyze all matrix tables at once.

Figure 7Five matrix tables presenting fitness metrics between synthetic model values and values obtained by the parameter estimation method using meteorological data from the year 2015 from the BEO site in Alaska. Matrix tables from left to right correspond to the normalized root-mean-squared errors for temperatures, liquid water saturations, and apparent resistivities as well as to the normalized Euclidian distances between synthetic (true) and estimated values of porosities and thermal conductivities. Each matrix value was normalized by dividing it by the matrix maximum value. The normalized values are indicated by tildes.

Download

To compare the match between all estimated and observational values within a single plot we calculated Euclidean norms for each case independently:

The index i indicates the case number (see Table 2). Then we applied k-means clustering analysis to identify groups of cases with a similar match between data and estimated parameters. We divided all cases into four classes shown in Fig. 8. Class I indicates the best cases that provide an accurate parameter estimation as well as accurate matches with the synthetic true measurements. Class II includes the cases that have accurate parameter estimates and less accurate matches with the measurements. Class III indicates all cases that have less accurate parameter estimates but accurate matches with the measurements. Finally, class IV includes the cases that showed the worst performance in terms of parameter estimates and the worst matches with the measurements. We summarized the results from Fig. 8 in Table 3.

Class I (see Table 3) suggests that sensors located below the active layer as well as increasing the number of sensors lead to more accurate parameter estimation. In contrast, case no. 4 (corresponding to the one ERT snapshot, 1ρ_a) suggests that already one ERT data snapshot could be enough for parameter estimation, while class II indicates that in general an increase in the numbers of monthly ERT snapshots is important for more accurate PE. However, increasing the number of monthly ERT snapshots leads to a less accurate match with measurements. These results are consistent with the results for simplified meteorological data with added noise (Fig. 5).

Class III includes six cases suggesting that if we have only soil moisture data available for PE, then we should expect less accurate soil property estimates. The last element in this class suggests that taking daily ERT snapshots improves the apparent resistivity match (Fig. 7, resistivity table) but does not improve {ϕ,k} estimates, where monthly ERT snapshots improve thermal conductivity convergence.

Class IV once again clearly indicates that measurements obtained below the active layer provide more accurate parameter estimates; however, they do not improve matches to measurements. This is mainly due to significant mismatch with ρ_a, which can be seen from the $\stackrel{\mathrm{̃}}{\mathrm{\Delta }{\mathit{\rho }}_{\text{a}}}$ matrix table in Fig. 8. At the actual site, the depth to the mineral soil can be deeper than 20 cm; not having sensors lower than 20 cm therefore limits the amount of data that can help to improve the convexity of the cost function in our case.

Table 3K-means analysis of the accuracy for each of the 13 cases.

Download Print Version | Download XLSX

Figure 8A k-means clustering analysis applied to the Euclidean norms of the normalized mean differences of estimated soil properties and the corresponding fit between calculated and observed values. Each color and marker represent a certain class as a result of the k-means clustering analysis.

Download

From Fig. 8 and Table 3 we know that the 6T case has the worst performance in terms of matching {ϕ,k}. Similar to the experiments with simplified meteorological data, the main difficulty for experiments with actual meteorological data is estimating thermal conductivity. The last matrix table ($\stackrel{\mathrm{̃}}{\mathit{\delta }k}$) in Fig. 7 shows that 6T6s_l8ρ_a(s) has the highest maximum and mean mismatch in thermal conductivity estimates. However, since ϕ estimates are a better match with their corresponding true values, the case 6T6s_l8ρ_a(s) falls into class III in Fig. 8, as opposed to case 6T, which falls into class IV. The highest mismatch in thermal conductivity values for the 6T6s_l8ρ_a(s) case suggests that collecting daily ERT snapshots improves the ρ_a match (Fig. 7, $\stackrel{\mathrm{̃}}{\mathrm{\Delta }{\mathit{\rho }}_{\text{a}}}$ matrix table) but does not improve estimated parameters, where monthly ERT snapshots improve thermal conductivity estimation.

To illustrate this, we plot values of estimated parameters and the corresponding response surfaces of the cost function for cases 10T10s_l8ρ_a and 6T6s_l8ρ_a(s) in Fig. 9. The PE method was able to match four out of five estimates almost perfectly and missed the k_p for the 10T10s_l8ρ_a case. The corresponding cost function has a visible minimum and clear convexity. In contrast to this, the 6T6s_l8ρ_a(s) case completely missed two estimates by converging on values outside the boundaries, and three other estimates do not converge to the desired cross section as well. The contour lines suggest that the corresponding response surfaces of the cost function do not have a well-defined global minimum.