## Abstract

Feedbacks between the land surface and the atmosphere, manifested as mass and energy fluxes, are strongly correlated with soil moisture, making soil moisture an important factor in land–atmosphere interactions. It is shown that a reduction of the uncertainty in subsurface properties such as hydraulic conductivity (*K*) propagates into the atmosphere, resulting in a reduction in uncertainty in land–atmosphere feedbacks that yields more accurate atmospheric predictions. Using the fully coupled groundwater-to-atmosphere model ParFlow-WRF, which couples the hydrologic model ParFlow with the Weather Research and Forecasting (WRF) atmospheric model, responses in land–atmosphere feedbacks and wind patterns due to subsurface heterogeneity are simulated. Ensembles are generated by varying the spatial location of subsurface properties while maintaining the global statistics and correlation structure. This approach is common to the hydrologic sciences but uncommon in atmospheric simulations where ensemble forecasts are commonly generated with perturbed initial conditions or multiple model parameterizations. It is clearly shown that different realizations of *K* produce variation in soil moisture, latent heat flux, and wind for both point and domain-averaged quantities. Using a single random field to represent a control case, varying amounts of *K* data are sampled and subsurface data are incorporated into conditional Monte Carlo ensembles to show that the difference between the ensemble mean prediction and the control saturation, latent heat flux, and wind speed are reduced significantly via conditioning of *K.* By reducing uncertainty associated with land–atmosphere feedback mechanisms, uncertainty is also reduced in both spatially distributed and domain-averaged wind speed magnitudes, thus improving the ability to make more accurate forecasts, which is important for many applications such as wind energy.

## 1. Introduction

The direct effects of subsurface heterogeneity have not been included in atmospheric studies to date. While the land surface and groundwater have historically been treated as simplified systems in atmospheric forecast and prediction models (Golaz et al. 2001; Kumar et al. 2006), early work by Chen and Avissar (1994), among others, has shown that soil moisture has a profound effect on local and mesoscale atmospheric processes. Work (e.g., Betts et al. 1996; Beljaars et al. 1996; Seuffert et al. 2002; Holt et al. 2006) has also shown that more advanced land surface model formulations and initialization, which generate more realistic soil moisture fields, result in better skill in mesoscale, regional, and local-scale weather forecasts. The reliance of these land surface parameterizations on accurate representations of surface soil moisture fields can be problematic because soil moisture is heterogeneous in space and time (Rubin and Or 1993; Or and Rubin 1993; Wendroth et al. 1999; Western et al. 2004; Famiglietti et al. 2008). Hydraulic conductivity, while highly variable and heterogeneous in space (several orders of magnitude), can be assumed to be static in time and has been shown to exhibit spatial correlation (Rubin 2003).

Many studies have connected soil moisture distribution and initialization to boundary layer atmospheric processes (e.g., Giorgi and Avissar 1997; Betts 2009). For example, seminal work by Chen and Avissar (1994) showed that soil moisture discontinuities, which can generate mesoscale atmospheric circulation patterns, have a strong influence on localized weather. The spatial distribution of precipitation is controlled by localized turbulent motion and mesoscale circulation, with the strongest rain intensities found at the mesoscale fronts. Later work (Patton et al. 2005) further explored a mechanism for atmospheric boundary layer (ABL) development over heterogeneous soil moisture fields in an idealized simulation in which they showed that total boundary layer turbulence kinetic energy increases significantly when the scale of heterogeneity is between four and nine times the depth of the ABL. The soil moisture–atmosphere connection has been confirmed through observational studies as well (e.g., Taylor et al. 2007). Recent studies (e.g., Chen and Hu 2004; Anyah et al. 2008; Jiang et al. 2009) show that land surface schemes that include groundwater dynamics further enhance forecast skill in shallow groundwater areas that would otherwise be moisture limited. Maxwell et al. (2007) and Maxwell et al. (2011) extend this to develop fully coupled three-dimensional subsurface flow and atmospheric models, and further show correlations between surface energy and moisture fluxes and the ABL.

Direct simulation of heterogeneous soil moisture conditions is difficult because soil moisture varies both temporally (e.g., Wendroth et al. 1999; Western et al. 2004) and spatially (e.g., Rubin and Or 1993; Or and Rubin 1993; Famiglietti et al. 2008) because of wetting and drying cycles and a dependence on topographic slope and hydraulic conductivity (*K*). Representing hydraulic conductivity in models presents one of the greatest challenges to the hydrologic and hydrogeologic community because it exhibits a high degree of uncertainty. Values of hydraulic conductivity can vary over several orders of magnitude within a single soil classification; for example, *K* values for a fine sand can fall within a range of 10^{−4}–10^{−7} m s^{−1} (Domenico and Schwartz 1998). Additionally, *K* is highly variable spatially, varying across several orders of magnitude over short distances. However, heterogeneity in hydraulic conductivity has been shown to be spatially correlated, which allows *K* to be represented statistically as a correlated random field (Rubin 2003). The uncertainty in the *K* correlated random field can be evaluated through multiple realizations in Monte Carlo ensemble simulations (Gelhar 1986) and reduced by assimilating observational data (Rehfeldt et al. 1992).

Ensemble, or stochastic, approaches are common in both the atmospheric and hydrologic–hydrogeologic sciences. However, the approaches used in each of the communities differ significantly. Atmospheric ensembles, from numerical weather prediction to climate change simulations, are commonly generated through perturbations of initial conditions and choice of model parameterization (e.g., Lewis 2005; Sutton et al. 2006; Leutbecher and Palmer 2008), while ensembles in hydrology and hydrogeology are motivated through uncertainty in input parameters such as rainfall distribution (Gabellani et al. 2007; Ridolfi et al. 2008) or spatial variability in hydraulic conductivity (e.g., Criminisi et al. 1997; Nowak et al. 2010). Additional work in hydrogeology has focused on reducing uncertainty in *K* through subsurface characterization (e.g., Vasco et al. 1997; Scheibe and Chien 2003). A common subsurface characterization approach in risk assessment, solute transport, and aquifer remediation studies employs Monte Carlo simulation ensembles to back-calculate *K* using observations of solute concentration or arrival times to condition realizations of the subsurface (e.g., Graham and McLaughlin 1989; Katul et al. 1993; James and Gorelick 1994; Harvey and Gorelick 1995; Yeh et al. 2005). Another approach directly conditions subsurface realizations to observed hydraulic conductivity data by assimilating point observations into a statistical representation of the subsurface (Maxwell et al. 1999; Rubin et al. 2010).

We apply a conditioning method whereby the distribution of *K* values in a correlated stochastic random field is controlled by enforcing “observed” point values drawn from a control simulation (CTRL) using a linear regression technique through which the stochastic random field honors both the observational data and the specified global statistics (Goovaerts 1997). Using a Monte Carlo simulation technique for both unconditioned and conditioned simulations, we show that *K*, saturation, latent heat flux, and wind speed magnitudes more closely honor hypothetical observed data and that improvements in atmospheric ensembles can be achieved by assimilating subsurface data.

## 2. Methods

We use the ParFlow hydrologic model coupled with the Weather Research and Forecasting Model (PF.WRF) to simulate subsurface, surface, and atmospheric conditions in a hypothetical 15 × 15 km basin. PF.WRF (Maxwell et al. 2011) couples the mesoscale WRF atmospheric model (Skamarock and Klemp 2008) with ParFlow, a three-dimensional variably saturated hydrologic model that simulates both subsurface and surface flow (Ashby and Falgout 1996; Jones and Woodward 2001; Kollet and Maxwell 2006), via mass and energy fluxes passed through the Noah land surface model (Chen and Dudhia 2001). Figure 1 shows a schematic cross section of the idealized model domain indicating which model component handles the subsurface, land surface, and atmosphere. ParFlow simulates flow in the subsurface by solving the Richards equation in three spatial dimensions with integrated surface flow routing using the kinematic wave approximation of the shallow water equations. The connection to evapotranspiration (ET) is simulated through moisture stress functions that limit the actual ET against potential ET. Additional details of ParFlow’s governing equations and the connections with the Noah land surface model are provided in the appendix.

To account for differences in model physics and parameterizations between a watershed model and an atmospheric model, the land surface–subsurface (ParFlow) and atmospheric (WRF) portions of the PF.WRF domain are set up separately. In the idealized simulations, horizontal discretization is the same for both the surface–subsurface and the atmospheric portions of the model domain—a 15 km × 15 km grid with 1000-m resolution. The subsurface portion of the domain was constructed with a vertical resolution of 0.5 m, a thickness of 5 m, and a 0.001 slope toward the west. The atmospheric portion of the domain is discretized in terms of 25 transient pressure levels defined fractionally from the land surface to the top of the atmosphere at 14 960 Pa, corresponding to a nominal vertical resolution of approximately 200 m near the surface and 800 m near the upper boundary of the domain at the top of the atmosphere—approximately 13 800 m above mean sea level.

The subsurface portion of the domain uses Neumann-type lateral and bottom boundary conditions and an overland flow boundary condition for the land surface, such that water enters the system only via rainfall and exits only via overland flow or evapotranspiration. Dynamic atmospheric forcing is passed from the atmospheric portion of the model in the form of mass and energy fluxes. Periodic lateral boundary conditions were imposed in all directions surrounding the atmospheric portion of the domain. The idealized model simulations are forced only by radiation inputs. The atmosphere is initialized with a slightly stable temperature profile given by *T*(K) = 300.0 − 0.005*z* where *z* is the height above the surface in meters, 50% relative humidity throughout the entire domain, and hydrostatic atmospheric pressure is based on a temperature of 300 K. Mean winds were not specified in the initialization, allowing winds within the domain to develop purely as a result of land–atmosphere feedbacks. We initialized soil moisture for these simulations by placing the water table at the bottom of the domain and applying 3 h of rainfall uniformly over the domain at a rate of 2 cm h^{−1}, using ParFlow in standalone mode (i.e., not coupled to WRF) before starting the fully coupled PF.WRF runs.

The simulation approach employed in this study is exploratory in nature, with a goal of understanding relationships between uncertainties in the subsurface and in atmospheric variables. This work was completed with a hypothetical test case, simulated actual conditions, and idealized initial and boundary atmospheric conditions at the field scale. Using a system with such simplifications allows us to explore the utility of this forecasting tool in a highly controlled setting in which we can isolate the effects of subsurface conditioning on the generation of land surface–controlled wind speeds. It has been shown in previous work (e.g., Famiglietti et al. 2008) that soil moisture varies over scales ranging from the field to the regional scale, and that field- and regional-scale heterogeneity in soil moisture affects mesoscale weather patterns (Chen and Avissar 1994; Patton et al. 2005). We ran this experiment using the Purdue–Lin microphysics parameterization (Lin et al. 1983) and repeated the entire experiment using the Morrison et al. (2005) two-moment microphysics parameterization scheme with substantially similar results.

Using PF.WRF, we ran four sets of Monte Carlo simulations. Each simulation comprised 10 realizations using different, yet statistically equivalent, heterogeneous *K* fields in the subsurface (please see the appendix for more details). Each simulation used stochastic random fields generated using the turning bands random field generator algorithm (Tompson et al. 1989) with an average (geometric mean) *K _{g}* of 0.01 m h

^{−1}, a variance of log

*K*of 2.25, horizontal correlation lengths of 5000 m, and vertical correlation lengths of 2 m for each realization. One additional realization, with its own random seed not included in the unconditional ensemble, was used to represent the CTRL conditions of the hypothetical domain. The conditioned sets of Monte Carlo simulations used the same statistical parameters and random seeds to generate the random fields as the unconditioned sets, but each was conditioned with an increasing number of data points (60, 120, and 200 points) drawn from the CTRL hydraulic conductivity field. The conditioning process, also described briefly in the appendix, requires that the values of hydraulic conductivity at sampled locations honor the CTRL value of

*K*. This results in a reduction in the variance around conditioning locations following Eq. (A13). Values of

*K*at conditioned locations in each ensemble member will be the same as the

*K*values at those locations in CTRL. Values of

*K*for cells in the vicinity of the conditioned locations in the ensemble members will be close to the

*K*values sampled from the CTRL case. Cells located more than one correlation length from any conditioning point have a variance equal to the field variance of 2.25. Parameterization selections and settings for these model runs are shown in Table 1.

We analyzed the model output data comparing the ensemble averages from the unconditioned (NO) case and each of the conditioned cases (CO60, CO120, and CO200) with the CTRL values for saturation, latent heat flux, and wind speed. These saturation and latent heat flux values are sampled from the surface soil layer. Wind speeds were calculated as cell-centered magnitudes from staggered-grid three-dimensional wind vector outputs. For the two-dimensional time slices shown in Figs. 2, 3, 5, and 7, we examine wind speeds in the pressure level closest to the surface. We evaluate uncertainty in the two-dimensional time slices at time *t* = 8 h using the mean squared residual *γ* between realizations and CTRL, calculated as

where *X* is the individual measurement for a realization and *n* is the number of realizations. This measure was used to quantify the residual between simulated and CTRL and to capture the variance within the ensemble in a single metric. We also plot domain-averaged time series for saturation, latent heat flux, and wind speed, shown in Figs. 4, 6, and 8. Root-mean-square (RMS) error between CTRL and the ensemble averages for each of the NO, CO60, CO120, and CO200 cases across the time series are shown in Table 2. Maximum *γ* and maximum variance within the time series also appear in Table 2.

## 3. Results and discussion

We focus our analysis on saturation, latent heat flux—variables that provide an indicator of surface conditions and energy fluxes as they relate to land–atmosphere feedbacks—and wind speed magnitude as the primary atmospheric variable of interest. We also focus on these variables as we expect the most direct and significant effect of conditional simulation on *K*, then saturation, then latent heat, and finally wind speed magnitude, spanning the subsurface to the atmosphere. Our analysis takes two parts: 1) we first examine *K*, saturation, and latent heat flux at the surface and wind speeds at the pressure level closest to the surface in two dimensions at a time slice of 8.0 h following the cessation of uniform rainfall, corresponding to the peak domain-averaged wind speed magnitude (as seen in Fig. 8; we examine the lowest elevation pressure level because, with a nominal vertical resolution of approximately 200 m, this is the level for which a wind forecast would be most relevant to wind energy applications); 2) we then analyze domain-averaged time series for saturation and latent heat flux averaged over the land surface and wind speed magnitude averaged over the entire atmospheric domain in the endpoint cases—NO and CO200. Results from the CO60 and CO120 cases fall within the bounds of these endpoints and are not shown.

### a. Hydraulic conductivity

The first variable we examine is hydraulic conductivity, which is directly affected by conditioning. The CTRL values for natural log of *K* are shown in Fig. 2a. As shown in Fig. 2b, *γ* values are highest in the NO case. In the CO60 case, values of *γ* decrease sharply in the area where *K* is conditioned (Fig. 2c). As expected, *γ* goes to zero at the conditioning points where the observed value of *K* is enforced. The spatial effect of the enforced *K* values [see Eq. (A12) in the appendix] can be clearly seen in the reductions of *γ* values in the vicinity of the conditioning points. With additional conditioning in the CO120 and CO200 cases (Figs. 2d,e), we see significant reduction in the *γ* values throughout the conditioned area. In the area to the east of the conditioning points, there is little to no improvement in *γ* for *K* as this area lies beyond the correlation length specified for the *K* stochastic random field. These results are illustrative of the method and are quite similar to what has been demonstrated in many other studies.

### b. Saturation

Saturation, shown in Fig. 3, is highly correlated with *K.* Similar to *K*, we observe a significant reduction in *γ* values between the NO case (Fig. 3b) and the CO60 case (Fig. 3c) and further reductions with increasing conditioning of the *K* random field in the CO120 and CO200 cases (Figs. 3d,e). The greatest reduction in *γ* occurs in the conditioned area, and little improvement is seen in the eastern portion of the domain (in the right side of each figure) where conditioning does not take place. Saturation is subject to lateral overland [Eq. (A3)] and three-dimensional subsurface flow [Eq. (A1)—i.e., infiltration], so as shown in Fig. 3, there are additional *γ* reductions seen over those seen for *K*.

The ensemble average saturation for the NO case overpredicts the CTRL in the domain-averaged time series (Fig. 4). In the CO200 case, saturation is slightly underpredicted by the ensemble average, but the RMS error is reduced by a factor of 4 (Table 2). The maximum variance and maximum *γ* value are also lower in the domain-averaged time series for the CO200 case.

### c. Latent heat flux

Latent heat flux, shown in Fig. 5, is heat transfer from the surface via evapotranspiration—a process that is limited by water availability and strongly correlated with saturation. As can be expected with this strong correlation, the behavior of latent heat flux closely resembles that of saturation and *K.* There are clear reductions of *γ* in the conditioned area with increasing conditioning of *K* (from the NO case through the CO200 case; Figs. 5b–e). However, because latent heat flux is a function not only of saturation but also of radiation flux, relative humidity, and wind (among others), its spatial correlation is not as restricted as saturation, which is more directly tied to the correlation lengths used to generate the *K* stochastic random field. Because of the averaging effect of atmospheric processes, we see small changes in *γ* values for latent heat flux in the eastern part of the domain outside the conditioned area, indicating that the effects of conditioning the subsurface may influence land–atmosphere feedbacks and weather patterns not only in the area where conditioning takes place but elsewhere as well.

Latent heat flux exhibits similar behavior to saturation in the domain-averaged time series (Fig. 6). The ensemble average overpredicts latent heat flux in the NO case but is very close to the CTRL values in the CO200 case. RMS error is reduced by a factor of 3 across the time series, and maximum variance and *γ* values decrease by factors of 3 and 4, respectively, between the NO and CO200 cases.

### d. Wind speed

Wind speed magnitude, shown in Fig. 7, exhibits the highest *γ* values where the wind speed is highest for the NO case. The strongest winds in the CTRL case (Fig. 7a) and in each ensemble average are predominantly from the west, defining clear up- and downwind directions. While the differences in *γ* values between the unconditioned and conditioned cases are not as dramatic on a domain-wide basis for wind as they are for land surface variables, the influence of subsurface conditioning is clear both in and out of the conditioned area, particularly in the high wind speed areas. Also noteworthy, especially in the CO200 case (Fig. 7e), is the reduction in *γ* values downwind of the strongest winds on the eastern part of the domain outside the conditioned area. The effects of subsurface conditioning are distributed across the domain and are not localized in the conditioned area, which is in contrast to the land-based variables. Indications of this effect are evident in the changes in *γ* values for latent heat flux with conditioning but are more clearly observed in wind speeds.

The domain-averaged time series for wind (Fig. 8) shows significant improvement in the ensemble averages between the NO and CO200 cases. The RMS error between the ensemble average and CTRL reduces by a factor of 2 from the NO to the CO200 cases. Contaminant transport studies employing conditioned Monte Carlo simulations (e.g., Maxwell et al. 1999) have used the ensemble variance as a measure of uncertainty, and improvements in the ensemble mean variance provide an indicator of reduction of uncertainty. The mean squared residual used in this work is a similar indicator, but in addition to providing a measure of the spread between ensemble members, it also provides a metric comparing ensemble member values to the CTRL values instead of the smoothed ensemble average. While the maximum variance between realizations and the ensemble average does not decrease appreciably with conditioning (in contrast to saturation and latent heat flux), the maximum mean squared residual does, indicating a greater likelihood that the ensemble members fall near the CTRL values. This is evident at the peak CTRL wind speed shown in Fig. 8 for the CO200 case where more than half of the realizations approach the peak at time 8.0 h—a significant improvement. Only one realization approaches the peak in the NO case.

## 4. Conclusions

Using a fully coupled subsurface-to-atmosphere model, we demonstrate that an atmospheric simulation ensemble can be generated with different realizations of subsurface *K*—a new finding. We further demonstrate that by conditioning *K* with an increasing number of observations, it is possible to reduce uncertainties in not only subsurface variables like saturation, but also in atmospheric variables such as wind; also a new finding. It has previously been established that ABL conditions are tightly coupled to soil moisture and latent heat flux from the land surface. It has also been previously established that soil moisture is a function of, among other variables, *K.* Through conditioning of *K* fields, we bridge these previous findings and demonstrate reduced uncertainty in the predicted soil moisture field, in latent heat flux, and in wind speed. The effects of conditioning the *K* field are evident in both spatially distributed cases and domain-averaged cases.

The reduction of uncertainty in *K* and the associated reduction of uncertainty in atmospheric variables are applicable to wind energy forecasts and to weather and atmospheric forecasts in general. Using a relatively small number of measurable observations, it is possible to reduce uncertainty in wind speed forecasts, which should prove useful in a wide range of atmospheric forecasting applications and climate change predictions.

The work we have presented is exploratory and intended to identify statistical relationships in the way uncertainty propagates between coupled fields in the subsurface and the atmosphere. To that end, we have completed a hypothetical test case with a small ensemble of semi-idealized model runs, and we have demonstrated these relationships. There are several facets of further work that could follow this study. These may include exploring the effects of coupled versus decoupled groundwater dynamics and perturbations in precipitation rates on the coupling of uncertainties in the subsurface and atmosphere to find the limits of the relationships we have shown here. Using a minimal number of realizations for our ensembles, we were able to show that uncertainty propagates through coupled fields in the subsurface, surface, and atmosphere. An expanded ensemble with many more realizations could improve the statistical representation of the uncertainties we have evaluated here and further improve forecast predictions. Finally, applying this technique to an actual site with nonidealized initial and boundary conditions could show whether or not this is a viable tool for operational weather forecasting.

## Acknowledgments

The authors gratefully acknowledge funding from the United States Department of Energy Wind and Hydropower Technology Program, Renewable Systems Interconnect Support. The authors also wish to acknowledge the peer reviewers of this article for their insightful comments.

### APPENDIX

#### Governing Equations and Statistical Techniques

While the equations solved by the ParFlow hydrologic model have been presented in detail elsewhere (e.g., Maxwell et al. 2011), a brief summary of the equations pertinent to the current work is presented here. The ParFlow hydrologic model accounts for pressure head in the subsurface and on the surface (Ashby and Falgout 1996; Jones and Woodward 2001; Kollet and Maxwell 2006). PF.WRF (Maxwell et al. 2011) accounts for moisture fluxes across the surface layer into the atmosphere using a formulation for evapotranspiration in the Noah land surface model. ParFlow calculates subsurface flow using the Richards equation in three dimensions according to Eq. (A1), with dependence on saturated hydraulic conductivity *K* shown in Eq. (A2) (Richards 1931):

where

In Eqs. (A1) and (A2), *S _{s}* is specific storage,

*S*is relative saturation,

_{w}*φ*is porosity,

*h*is pressure head,

*q*is a source–sink term that includes precipitation and evapotranspiration,

_{r}**K**

*(*

_{s}**x**) is the saturated hydraulic conductivity tensor,

*k*is the relative permeability,

_{r}**v**is the subsurface flow velocity, and

**q**is the specific volumetric flux.

ParFlow represents overland flow by the two-dimensional kinematic wave equation [Eq. (A3)], which also depends on *K.* The kinematic wave formulation is part of the overland flow boundary condition (Kollet and Maxwell 2006)

where **v**^{sw} is the depth-averaged surface water velocity in two dimensions, *h* is the surface ponding depth and pressure head at the surface, and **k** is the vertical unit vector. The greater of pressure head *h* and zero is represented by ‖*h*, 0‖ in Eq. (A3). Equation (A3) connects flow processes in the subsurface with those at the surface.

Evapotranspiration is calculated in the Noah land surface model as a function of vegetation fraction and potential evaporation as follows:

where *E* is the rate of evapotranspiration, fx is an empirical coefficient, *f*_{veg} is the vegetation fraction, and *E*_{pot} is potential evaporation, calculated based on atmospheric conditions from the WRF boundary layer parameterization (Ek et al. 2003). The atmospheric experiment shown in this work is conducted over bare soil, so vegetation fraction is zero. Thus, Eq. (A4) reduces to

The quantity *F* is parameterized as follows:

where *S*_{res} is residual saturation. Relationships between pressure head and relative saturation and relative permeability are parameterized using the van Genuchten functions (van Genuchten 1980).

Each of Eqs. (A1) through (A6) depend on *K,* either directly or through functions of saturation and pressure head. The parameterization of *K* is crucial for calculations of movement of water in the subsurface and on the surface and for calculations of moisture and energy fluxes from the land surface to the atmosphere. Because of the high level of uncertainty associated with hydraulic conductivity, it is not possible to generate a representation of the subsurface that is both complete and accurate. Instead, we generate a statistical representation of subsurface hydraulic conductivity as a spatially correlated stochastic random field, which can be conditioned through subsurface characterization and assimilation of observed data. We use a natural log transform of hydraulic conductivity where

Perturbations in the hydraulic conductivity field take the form

where angle brackets represent an expected or mean value, and *y*′ is the magnitude of the perturbation. The geometric mean of hydraulic conductivity *K _{g}* is taken to represent the expected value of

*Y*:

The expected value of the perturbation is zero:

however, the square of the expected value of the perturbation is the variance:

(Rubin 2003).

The anisotropic correlation structure of hydraulic conductivity values between any two points is defined by the covariance

where *r* is the separation distance, *λ* is the correlation length, and the subscript *d* represents the spatial direction (Rubin 2003). The correlated random field is generated using the turning bands random field generator algorithm (Tompson et al. 1989). The random field generator is initialized with a random number seed. Multiple realizations within an ensemble, each with a different random seed, will have different, but statistically equivalent, random fields of hydraulic conductivity.

We condition the hydraulic conductivity fields by assimilating hydraulic conductivity values sampled from a control or “truth” realization that is independent of the model ensembles. Conditioning is accomplished by enforcing sampled values from the control simulation at their corresponding locations within the stochastically generated correlated random field. Using simple kriging, the perturbation in the hydraulic conductivity field, *y*′ in Eq. (A8), is tied to observed data values via a weighted correlation structure. At each point in the correlated random field, for *n* observed data points,

The summation in Eq. (A13) relates the perturbation, *y*′ in Eq. (A8) to observed data. The weight, *ξ*, is determined by calculating the minimum error variance at each point in the correlated random field for *n* observed data points:

where is the field variance [Eq. (A11)] and *C* is the covariance as defined in Eq. (A12) (Goovaerts 1997; Rubin 2003). The simple kriging system of equations incorporates the global statistics defined for the random field generator, and it relates every point in the correlated random field to all observed data points in calculating the weight factor based on covariances. In this way, the correlated random field representing hydraulic conductivity honors assimilated data point values while still also honoring the specified global statistics in each conditioned realization.

## REFERENCES

## Footnotes

^{*}

Additional affiliation: Integrated Groundwater Modeling Center, Colorado School of Mines, Golden, Colorado.