Sequential and joint hydrogeophysical inversion using a field-scale groundwater model with ERT and TDEM data

. Increasingly, ground-based and airborne geophysical data sets are used to inform groundwater models. Re-cent research focuses on establishing coupling relationships between geophysical and groundwater parameters. To fully exploit such information, this paper presents and compares different hydrogeophysical inversion approaches to inform a ﬁeld-scale groundwater model with time domain electromagnetic (TDEM) and electrical resistivity tomography (ERT) data. In a sequential hydrogeophysical inversion (SHI) a groundwater model is calibrated with geophysical data by coupling groundwater model parameters with the inverted geophysical models. We subsequently compare the SHI with a joint hydrogeophysical inversion (JHI). In the JHI, a geophysical model is


Introduction
Over the last decade, interest in geophysical methods for hydrogeological site characterization has been increasing (Vereecken et al., 2004;Hubbard and Rubin, 2000). This is due to the ability of geophysical methods to provide models of subsurface properties with a high spatial resolution, which are difficult to obtain from sparse borehole information. Worldwide, significant resources are spent on the collection of regional geophysical data sets. Examples include airborne electromagnetic (AEM) surveys in Denmark, covering nearly 60 % of the country for mapping the spatial extent and assessing the vulnerability of aquifers (Thomsen et al., 2004), and AEM surveys to map saltwater intrusion in the USA, Australia, Germany and the Netherlands (Langevin et al., 2003;Fitzpatrick et al., 2009;Faneca Sànchez et al., 2012;Burschil et al., 2012). In addition, smaller-scale surveys have been conducted using a variety of geophysical techniques such as ERT (electrical resistivity tomography, Kemna et al., 2002), induced polarization (Slater, 2007) and magnetic resonance sounding (Legchenko and Valla, 2002).
Published by Copernicus Publications on behalf of the European Geosciences Union. 1 2 Fig. 1 SHI (left), CHI (middle) and JHI approach (right). π and γ respectively indicate the geophysical and groundwater model parameters, where the 3 arrows represent parameter updating until a minimum data and/or constraint misfit is achieved. π: geophysical model parameters; γ: groundwater model 4 parameters. 5 Fig. 1. SHI (left), CHI (middle) and JHI approach (right). π and γ respectively indicate the geophysical and groundwater model parameters, where the arrows represent parameter updating until a minimum data and/or constraint misfit is achieved. π: geophysical model parameters; γ : groundwater model parameters.
A major challenge is to fully exploit the information content of geophysical data sets, as geophysical techniques do not measure hydrological subsurface properties directly. A geophysical inversion and petrophysical relationships are needed to estimate hydrogeological parameters and state variables from the geophysical data sets. For this reason, the inclusion of geophysical data into a groundwater model is not straightforward. Previous studies have used different approaches to inform groundwater models with geophysical data.

Hydrogeophysical inversion approaches
Hydrogeophysical inversion approaches can be subdivided into sequential hydrogeophysical inversion (SHI), coupled hydrogeophysical inversion (CHI) and joint hydrogeophysical inversion (JHI) (Hinnell et al., 2010). The workflow associated with these 3 methods is shown in Fig. 1. 1. In a SHI, geophysical data is separately inverted to estimate the distribution of a geophysical property (e.g., maps of electrical resistivity). Estimated geophysical property maps are subsequently used to derive the structure of the subsurface or to estimate dynamic state variables such as solute concentrations and water content. For the latter, petrophysical relationships (Archie, 1942;Topp et al., 1980) are needed to convert a geophysical property to a hydrological state variable.
Note Fig. 1 only shows an SHI in which inverted geophysical parameters are coupled with the static input structure of a hydrological model; SHI by coupling dynamic state variables would typically require coupling inverted geophysical parameters with hydrologic model output.
2. In a CHI, simulated state variables of a hydrological model are transformed to a geophysical parameter distribution using a petrophysical relationship. Subsequently, geophysical forward responses are simulated that can be compared with collected geophysical observations. In this approach, the geophysical inversion process is coupled with the hydrological model and a single objective function is minimized that comprises both a geophysical and a hydrological component.
3. In a JHI, a simultaneous inversion for multiple geophysical and/or hydrological models is undertaken to exploit differences in parameter resolution for different data sets. In the JHI discussed in this paper, input parameters of a hydrological and geophysical model are simultaneously estimated using parameter coupling constraints to account for observed correlations between geophysical and groundwater model parameters (e.g., petrophysical relationships).
Examples of SHI applications include the use of geoelectrical methods, electromagnetic methods and ground penetrating radar (GPR) to monitor changes in soil water content or solute concentrations with time (Binley et al., 2001;Cassiani et al., 2006;Day-Lewis et al., 2003;Huisman et al., 2003;Kemna et al., 2002;Knight, 2001;Looms et al., 2008). Of particular interest is the SHI framework presented by Dam and Christensen (2003) in which inverted electrical resistivities are used to estimate hydraulic conductivity fields of a groundwater model. As will be explained later, our JHI approach shows many similarities with this framework. Examples of CHI applications include the estimation of vadose zone parameters with electrical resistivity and GPR measurements (Hinnell et al., 2010;Kowalsky et al., 2005;Lambot et al., 2006Lambot et al., , 2009, the estimation of hydraulic conductivity fields with electrical resistivity data (Pollock and Cirpka, 2012) and the estimation of soil properties with relative gravimetry and magnetic resonance sounding data (Christiansen et al., 2011;Herckenrath et al., 2012a). These studies cover a relatively small spatial scale compared to field-scale groundwater models. Applications of a CHI on a more regional scale can be found in (Bauer-Gottwein et al., 2010;Herckenrath et al., 2012b). JHI methods have been developed to use multiple geophysical methods for estimating soil properties (Vozoff and Jupp, 1975;Linde et al., 2006a;Behroozmand et al., 2012) (Hyndman and Gorelick, 1996;Chen et al., 2006;Linde et al., 2006b;Herckenrath et al., 2012a). In comparison with a SHI, the main strength of a CHI is to use a hydrological model to provide an inversion framework for the geophysical data and constrain the geophysical inversion process. This brings two main advantages, as (1) measurement errors and parameter uncertainties associated with the independent geophysical inversion are not propagated directly to the hydrological model and (2) no extensive regularization (e.g., smoothness constraints) is needed to stabilize the geophysical inversion process (Menke, 1984). These regularization constraints do not necessarily reflect the hydrological conditions (Day-Lewis et al., 2005;Chen et al., 2006;Slater, 2007) and in a CHI these are substituted by spatial correlation structures provided by a hydrological model. In simple words, the hydrological model provides an interpretation framework for the geophysical data.
The advantage of joint inversion/JHI with respect to SHI is the exploitation of parameter resolution differences for different data types (Linde et al., 2006a). Concerns pertaining to joint inversion for multiple geophysical methods are mainly related to observation weighting strategies and the transfer of correlated measurement error. When different model types and setups are used, as in a JHI, transfer of conceptual model errors is an additional problem.
In this study we confront the SHI and JHI approaches, without focusing on the JHI and CHI comparison, as these latter methods cannot easily act as substitute for one another. This is due to the different nature of coupling between the hydrologic and geophysical model. In a CHI, hydrological simulations are coupled with geophysical models while JHI couples input parameters.
However, JHI and CHI share similar concerns regarding the propagation of hydrological conceptual model errors, the definition of petrophysical relationships and the assignment of weights for various data types. The propagation of hydrological conceptual errors to the geophysical model differs for SHI, CHI and JHI. In a SHI no conceptual hydrological errors propagates into the geophysical model, but this will be the case in CHI and JHI. The difference between the latter two methods is that CHI generally employs a single conceptualization for both the geophysical and hydrological model while JHI allows different conceptualizations for the geophysical and hydrological model. Further discussion of this topic is beyond the scope of this paper, as this problem will be highly dependent on considered models, available data sets, the purpose of the hydrogeophysical inversion process and the employed petrophysical relationships.

Petrophysical and geometric relationships
Any hydrogeophysical modeling approach (SHI, CHI or JHI) depends on coupling geophysical and hydrological models by implementing coupling relationships between geophysical parameters with hydrological model parameters or hydrological model simulations. Such coupling relationships can be sub-divided in different groups. In this paper, we consider petrophysical and geometric relationships. Well-known petrophysical relationships are Archie's law (Archie, 1942) and the Topp equation (Topp et al., 1980), that respectively link electrical resistivity and relative electrical permittivity with hydrological properties such as porosity and water content. In the context of field-scale groundwater modeling, relationships between hydraulic conductivity and geophysical properties would be of particular interest. Research shows that such relationships exist, including a log-linear correlation between hydraulic conductivity and electrical resistivity (Purvance and Andricevic, 2000;Niwas and de Lima, 2003), the dependence of chargeability on clay-content (Slater, 2007) and the estimation of hydraulic conductivity from magnetic resonance sounding data (Vouillamoz et al., 2008). Typically, petrophysical relationships are site-specific and are established based on field observations. Site-specific relationships might be extrapolated for hydrogeological units within the same sedimentary basin, as previous studies showed the importance of taking into account geological properties for obtaining a petrophysical relationship (Prasad, 2003;Slater, 2007).
Geometric relationships comprise the use of structures derived from geophysical models to identify spatial geological information used in hydro(geo)logical models. An example is given in Burschil et al. (2012), in which AEM, seismic reflection and borehole data is used to define the hydrostratigraphy of a groundwater model for a complex glacially affected island. Hydrostratigraphy can be estimated as part of hydrogeological model calibration (Passadore et al., 2012), in which geometric constraints can be used to tie the hydrostratigraphy of a groundwater model with a geophysical model. This can be relevant for the definition of confining units and saltwater intrusion models, where aquifer thickness and bathymetry are important properties (Carrera et al., 2010).

Aim of this study
Hydrogeophysical inversions are generally used for smallscale studies. Given the developments in geophysical data collection for regional groundwater exploration and available work on petrophysical relationships, we aim to extend the use of hydrogeophysical inversions for field-scale and regional groundwater models. For this purpose, we implement and compare JHI and SHI for a field-scale groundwater model with TDEM and ERT data. The study faces a number of specific challenges: 1. The conceptual framework of groundwater models is prone to large uncertainties (Refsgaard et al., 2006), due to their scale, limited data availability and the use of many simplifying assumptions associated with the geological model and boundary conditions.
2. The sub-surface volumes represented by groundwater and geophysical models can be very different, which limits using a single conceptualization for both models.
3. Some subsurface processes and/or compartments are included in the geophysical or hydrogeological model only and are not represented in the other model.
4. The accuracy of the relationship between geophysical and groundwater model parameters is difficult to determine.
5. Computational burden and large number of estimated parameters.
Based on the first three issues, geophysical and hydrogeological models usually require different conceptualizations to achieve an acceptable data fit. This flexibility cannot be incorporated when the geophysical model is completely constructed from hydrological model in-or output as in many CHI studies. The strength of coupling between the geophysical and groundwater models is difficult to determine and can be based on the assumed accuracy of the (petro)physical relationships between geophysical and groundwater properties. This accuracy can be estimated from correlating geophysical models with available groundwater data (e.g., pumping tests, borehole data, and lab tests). In a SHI the strength of coupling constraints can be either based on geophysical parameter resolution or the accuracy of the petrophysical relationship.
The fifth challenge is related to the large computational burden associated with groundwater models and inversion of geophysical models. Due to the computational burden, parameter estimation is typically performed using local, gradient-search algorithms (Doherty, 2010) instead of global search algorithms like Markov-Chain Monte Carlo based methods (Vrugt et al., 2009). Gradient-search algorithms, such as the Levenberg-Marquardt method, do not always find the true global minimum of the objective function surface due to multiple local minima in parameter space, discontinuous first derivatives, curved multidimensional ridges and parameter surrogacy (Vrugt et al., 2008). Initial parameter values are therefore extremely important when using local, gradient-search techniques.
The final challenge refers to correlated geophysical measurement errors that can be caused by existing infrastructure (e.g., power lines, buried pipes), neglecting 3-D effects in the geophysical model (Bauer-Gottwein et al., 2010) and the application of inaccurate/limited instrument filters when processing geophysical data (Efferso et al., 1999). Characteristics of correlated noise are location-specific and different for the various types of geophysical methods and therefore difficult to quantify. We do not consider correlated measurement error in this paper. An example of how correlated measurement error propagates in a CHI is provided in Hinnell et al. (2010) and Herckenrath et al. (2012a).
To meet the previous mentioned challenges, we implement a SHI and JHI in which geophysical model parameters are tied to groundwater model parameters by adding parameter coupling constraints. These parameter coupling constraints can be imposed to subsets of parameters to ensure enough flexibility to fit the different types of observation data, while the imposed strength of the parameter coupling constraints reflects the quality of the relationship between model parameters that can be derived from field data or geophysical parameter resolution. Finally, these parameter coupling constraints are compatible with standard inversion methods used for groundwater and geophysical models. The presented SHI-approach is similar to Dam and Christensen (2003), whereas the JHI is similar to an inversion methodology used by Doherty and Johnston (2003), which differs from standard joint inversion approaches, as input parameters are not shared by multiple models but coupled through additional regularization constraints.
Section 2 provides a theoretical background for the applied SHI and JHI. Section 3 shows the application of both the SHI and JHI for a synthetic groundwater model with time domain electromagnetic (TDEM) data. The implementation of a JHI and SHI for a real-world groundwater model and geoelectric data (ERT) is described in Sect. 4. Results are given in terms of parameter estimates, parameter error, model misfit and computational burden. The paper concludes with a summary of the benefits, disadvantages and limitations associated with the presented coupling procedures.

Methodology
This section provides a mathematical summary of a SHI and JHI.

Sequential hydrogeophysical inversion (SHI)
The SHI starts with a geophysical inversion. Consider a data set of geophysical observations assembled in vector d g : The symbol ρ denotes the geophysical observations, e.g., apparent resistivities. Subscript N g is the number of available geophysical observations and e g denotes the geophysical measurement error. The geophysical model parameters that are estimated are assembled in vector π : π = (r 1 , ., r M r , t 1 ., t M t ) T . (2) In this paper π contains a number of layer thicknesses (t) and layer resistivities (r) for a 1-D electrical resistivity The SHI starts with a geophysical inversion in which geophysical parameters in π are estimated by fitting the geophysical observations in d g . For this purpose we follow a well-established iterative least-squares inversion approach (Tarantola and Valette, 1982;Menke, 1984).
According to Auken and Christiansen (2004), the inversion problem can be written as ( In the geophysical inversion, a geophysical forward model is used to calculate apparent resistivities for the electrical resistivity model defined in π. G g is the Jacobian comprising the partial derivatives of d g with respect to the geophysical parameters in π. Furthermore, four types of regularization constraints are used in the inversion: prior parameter constraints, prior depth constraints, vertical constraints and lateral constraints. These result in four additional operators I, P h , R p and R h and contribute to the total geophysical observation error e g . The implementation and derivation of these constraints is explained in detail in Auken and Christiansen (2004). δπ prior , δπ h-prior , δr p and δr h express the deviation with respect to the expected value for the prior parameter constraints, prior depth constraints, vertical constraints and lateral constraints. e prior , e h-prior , e p , and e h are the errors associated with these constraints. More compact Eq. (3) is In the geophysical inversion the following objective function is minimized by updating π, (5) where ϕ prior , ϕ h-prior , ϕ Rp , and ϕ Rh represent the objective function component for the additional parameter constraints as defined in Auken and Christiansen (2004). The posterior standard deviation of the estimated geophysical parameters is calculated based on a post-calibrated parameter covariance matrix, defined as where C g defines the parameter covariance matrix. Posterior parameter standard deviations are subsequently calculated as the square root of the diagonal elements of C gest using STD(π est ) = C gest (s, s), where π est represents the final geophysical parameter estimate and s = 1, 2, . . . , M g . Next, we consider a set of groundwater observations that are listed in vector d h , subscript N h indicates the number of groundwater observations represented by h, which can include head data and observed water fluxes. e h defines the measurement errors on the groundwater data. The groundwater model parameters are listed in vector where M h represents the number of groundwater parameters; in this paper these parameters represent hydraulic conductivities and thicknesses of geological layers. An iterative least squares approach is used to estimate the parameters listed in γ . For the groundwater data we write where G h is the Jacobian containing all partial derivatives associated with the groundwater forward mapping. The second step of the SHI is to calibrate the groundwater model using the traditional data in vector d h and a number of estimated geophysical model parameters π est together with their posterior standard deviations. When a petrophysical relationship is used, π est is first transformed to another property (e.g., hydraulic conductivity). This yields an additional set of hydrogeological observations comprised by vector s h , where N s is the number of transformed geophysical parameters, p, that are used as additional observations to constrain the groundwater model parameters. These observations are connected to the groundwater model parameters as given in Eq. (12): where P s is a matrix with the dimensions of γ and N s , containing ones for the groundwater model parameters that are constrained by the estimated geophysical parameters in s h and zeros for the groundwater model parameters that are not constrained. e s represents the posterior standard deviations associated with the geophysical parameters. This approach is analogous to the use of the prior parameter constraints in the geophysical inversion. The hydrogeological inverse problem can therefore be described as or more compact as with parameter update where C h ' is the joint observation error comprising the error covariance matrix C h for the hydrogeological observations and C s for the geophysical observations. Equation (15) minimizes the objective function ϕ SHI defined as Parameter uncertainty is calculated using a posterior parameter covariance matrix as described by Eq. (7). Note the SHI is equivalent to the method described in Dam and Christensen (2003), except for the definition of e s .

Joint hydrogeophysical inversion (JHI)
In a SHI the strength of coupling between the geophysical and groundwater model is based on e s , which in our implementation depends on geophysical parameter resolution only. Another coupling strategy would be to define the strength of coupling based on the accuracy of established petrophysical relationships.
In contrast to the SHI, JHI performs one single inversion for both the geophysical and the hydrogeological model. For this purpose, the parameters of both models are assembled in vector m, m = (π 1 , π 2 , ., π M g , γ 1 , γ 2 , ., γ M h ) T .
We introduce a number of coupling constraints between the geophysical and hydrogeological parameters that are connected to the true model as where e c denotes the error associated with the coupling constraint. Because the coupling constraints link different estimated parameters, e c is unknown and has to be defined by the user. Its definition depends upon the assumed error of the coupling constraint. e c plays a key role in the JHI framework and its value can be estimated from available field data that was used to establish a relationship between a groundwater and geophysical parameter. In Slater (2007) correlation plots are provided between geophysical properties and hydraulic properties. The correlation measure of such analyses can be used to estimate e c .
Operator P c can have many forms. For example, if we introduce two coupling constraints that set the groundwater model parameters γ 1 and γ 2 (geological layer thicknesses), equal to respectively π 1 and π 2 (e.g., geophysical model layer thicknesses), Eq. (18) takes the following form: Note that for petrophysical relationships between π and γ , δr c in Eq. (18) often has a nonzero value. An example will be provided in the case study section. Coupling constraints between π and γ need to be linear for the current implementation of the JHI.
Combining Eqs. (4) and (10) with the coupling constraints in Eq. (18), we obtain the formulation for the JHI: which can be written more compactly as Many of the entries in Jacobian G are equal to 0 as some of the hydrogeological parameter estimates are not affected by the geophysical observation and constraints and vice versa. The joint observation error e is denoted by covariance matrix C : The model estimate becomes which minimizes the objective function where φ h is the hydrogeological data misfit, ϕ g the geophysical data misfit and φ c the objective function term associated with the coupling constraints. φ c acts as an additional regularization term mutually constraining the geophysical and groundwater parameters. A similar approach can be found in Doherty and Johnston (2003), who estimate parameters of multiple watershed models.

Implementation
The SHI and JHI are applied for two cases. The first case combines a synthetic groundwater model and a synthetic TDEM data set. The second case combines a real-world groundwater model and a field ERT data set.
To generate the geophysical forward responses for the TDEM and ERT, EM1DINV (HGG, 2008) is used. EM1DINV is also used to generate a forward response for the ERT data (Auken and Christiansen, 2004). The geophysical model that is estimated for the TDEM is a 1-D resistivity model (Fig. 2b), in which typically a number of layer thicknesses and layer resistivities are estimated. For the ERT data, neighboring 1-D resistivity models (Fig. 8a) are tied together by lateral constraints (Auken and Christiansen, 2004).
The groundwater model in the synthetic example is implemented in Matlab (PDE-tool). For the real-world model MODFLOW (Harbaugh et al., 2000) is used. More details about the groundwater models and geophysical data are given in the next section.

Setup
The first application of the JHI and SHI considers a synthetic cross-sectional groundwater model and a TDEM sounding. As part of the geophysical inversion a TDEM forward model is used. This forward model is based on Ward and Hohmann (1988) and includes the modeling of low-pass filters (Efferso et al., 1999) and the turn-on and turn-off ramps described in Fitterman and Anderson (1987).
The groundwater model in the synthetic example consists of two layers, similar to the geological setup of the field study we discuss in Sect. 4. The upper layer, with a thickness D clay , is considered to be clayey sand with hydraulic conductivity K clay [m s −1 ]. The second layer represents limestone with hydraulic conductivity K lime . Different values are generated for these properties as will be explained below. Constant heads are applied as boundary conditions (right: 1 m; left: 0 m); in the middle of the model domain a river is assumed to be located with a fixed head of 0. This results in flow from left to right and flow towards the river. From this realization we extract a number of groundwater observations, comprising 4 head and 2 flux measurements that are shown in Fig. 2a. The groundwater parameters (γ ) that need to be estimated include the hydraulic conductivity of the limestone (K lime ) and the clay (K clay ) and the thickness of the clay (D clay ). Due to the parameter cross-correlation between K clay and D clay , an additional flux measurement for the limestone is included, which is not available for most real-world modeling studies. Typically D clay is not estimated when calibrating a groundwater model, due to its correlation with K clay . This parameter was chosen to illustrate the use of a JHI and SHI, in which the hydrostratigraphy of a groundwater model is coupled with a geophysical model.
For the synthetic study we assume the availability of one TDEM sounding. The parameters of the geophysical model (π ) that are estimated comprise one layer thickness (t 1 ) and electrical resistivities for layer 1 and 2 (r 1 and r 2 ) using 30 synthetic apparent resistivity observations. The simplified 1-D description of the geophysical model is used because of the negligible effect of the water table variation and unsaturated zone thickness in the model, compared to the geometry of the model and the TDEM resolution.
In summary, 6 parameters are estimated, 3 for the geophysical model and 3 parameters for the groundwater model. To test the SHI and JHI, we generate 250 observation realizations of hydrogeological data (heads and fluxes) and geophysical data (apparent resistivities) by adding uncorrelated measurement error to a model-generated truth. For every realization different values for K clay , D clay , r 1 and t 1 are generated, each representing a model generated truth. The generation of log10 K clay [m s −1 ] and D clay [m] values employed mean values of respectively −5 and 25 m with a standard deviation of respectively 0.1 and 0.1 m. Subsequently values of r 1 and t 1 are generated based on the equations in the second column of Table 2, including a random component with a standard deviation, e corr , that defines the level of correlation between the geophysical and groundwater model parameters. Measurement error is then added to the simulation results of each parameter realization, employing a standard deviation (e h ) of ±2 cm for the head observations and ±30 % for the flux measurements. The measurement errors added to the TDEM data have a standard deviation (e g ) of ca. ± 3 % of the measurement value and are based on a real-world TDEM sounding.
The TDEM measurement error does not only reflect the standard deviation of the data stack and includes an additional error component to take into account 3-D effects and imperfect instrument specifications (e.g., filters, wave form of the applied pulses). This additional error component will typically yield correlated measurement errors. For example, Efferso et al. (1999) provide the effect of different low pass filters on the TDEM forward response. In this research, however, we do not investigate correlated errors and thus add uncorrelated measurement error to the TDEM data to be consistent with the Gaussian assumptions of least-squares inversion theory (Tarantola, 2005). Different starting parameters are used for the calibration of the geophysical and groundwater model with each observation realization.

Geometric and petrophysical relationship
To perform the JHI and SHI two types of constraints are employed between the groundwater and TDEM model, a geometric and a petrophysical constraint. Both relationships are defined in Table 2. The geometric constraint applies to the depth of the clay layer (D clay ) and the thickness of the first layer in the TDEM model (t 1 ).
The petrophysical coupling constraint applies to the hydraulic conductivity of the upper layer of the groundwater model (K clay ) and the electrical resistivity of the first layer in the TDEM model (r 1 ). This constraint applies a relationship between the logarithmic values of hydraulic conductivity and electrical resistivity (Niwas and de Lima, 2003;Slater, 2007). The petrophysical relationship in Table 2 was arbitrarily chosen, but implies a decreasing hydraulic conductivity for a decreasing electrical resistivity, as hydraulic conductivity and electrical resistivity decrease for increasing clay content. A typical hydraulic conductivity for clay is 10 −5 m s −1 (Fetter, 1994) and 10 1 m is a representative electrical resistivity (Kirsch, 2006), which results in an expected value of −6 for the petrophysical coupling constraint. Note that this is an extremely simplified relationship between hydraulic conductivity and electrical resistivity.
In a first configuration of the synthetic study, we generate realizations of "true" parameters, using a standard deviation (e corr ) of 0.01 for the petrophysical relationship and a standard deviation 0.05 (e corr ) for the geometric relationship. In a second configuration, we apply larger e corr values of respectively 0.1 and 0.1. As the parameter coupling in the SHI can be very strong for well-resolved geophysical parameters, this second configuration is used to test whether or not the SHI results in worse groundwater parameter estimates when correlation between groundwater and geophysical parameters is relatively weak.

SHI
The SHI starts with a geophysical inversion for the TDEM data after which the estimated electrical resistivity model, π est , is used as an observation in the calibration process of the groundwater model. In this case π est comprises the estimated values for t 1 and r 1 , which we employ to constrain the groundwater model parameters D clay and K clay . For the weights of these constraints (Dam and Christensen, 2003) Hydrol. Earth Syst. Sci., 17, 4043-4060, 2013 www.hydrol-earth-syst-sci.net/17/4043/2013/   Table 2. recommend e s values of 10 −2 -10 −1 for coupling hydraulic conductivities and well-resolved electrical resistivities and values of 10 1 -10 2 for poorly resolved electrical resistivities. We employ values based on the posterior standard deviation of the geophysical parameters, obtained with the geophysical inversion, to honor the resolution level of parameters inferred from geophysical data and constraints.

JHI
For the JHI we use the same type of coupling constraints for the same geophysical and hydrological parameters. However, now the geophysical parameters are also part of the inversion and Eq. (18) is used for the coupling constraints. For this application Eq. (18) becomes where the expected value for the geometric constraint between D clay and t 1 is 0, whereas the petrophysical relationship between log10(K clay ) and log10(r 1 ) is 6. The JHI is undertaken for varying values of e c , as defined by the values in Table 2. This range is comparable with the recommended range for e s in Dam and Christensen (2003). The value of e c reflects the strength of the coupling relationship. An e c of 0.01 means the assumed error of the coupling relationship has a standard deviation of 0.01, marking a strong coupling relationship compared to an implementation employing and e c of e.g., 10. For the synthetic study the weight associated with the coupling constraints is varied, by changing this standard deviation. Table 2   were chosen to cover a JHI with weak coupling constraints and a JHI assuming e c values of similar magnitude compared to the standard deviations, e corr , that were used for generating the correlated "true" parameters.

Results
First a JHI is conducted for the groundwater and the geophysical model. This was done using 250 observation realizations and different parameter starting values. 7 JHI simulations are performed using an increasing strength of coupling between the TDEM and groundwater model (Runs 1-7). To generate correlated "true" geophysical and groundwater model parameters, standard deviations e corr of 0.01 and 0.05 are respectively used for the petrophysical and geometric constraint. Run 1 represents a JHI with a very small weight (i.e., large e c ) for the coupling constraints representing an independent inversion in which the groundwater model is not informed with the TDEM model and vice versa. Figure 3 shows all the parameter estimates pertaining to the JHI Run 1-7 for 250 realizations, expressing how well parameter estimates compare with the "true" parameter values that were generated. Parameter errors in Fig. 3 are given as a percentage with respect to the "true" parameter value. For JHI Run 1 parameter errors are up to 100 % for K clay and K lime and up to 40 % for D clay .
Geophysical parameter r 1 is well-resolved and shows errors of less than 7 %. t 1 and r 2 show errors of respectively 40 and 200 %.
The strength of the coupling constraints is subsequently increased using smaller values for e c (Table 2) in JHI Runs 2-7. The blue dashed lines in Fig. 3 shows how parameter estimates react as a result of the stronger coupling constraints. A large and rapid reduction of error can be observed for K clay showing an error decrease from 100 % to about 10 %. Estimates for D clay do not improve and errors remain at a value of up to about 40 %. Geophysical parameter errors are fairly constant for Runs 1-7, except for a slightly increasing number of realizations showing larger errors for parameter r 1 and t 1 in JHI Runs 6 and 7 in which the coupling constraints have the largest weight. Figure 4 provides the data fit for the different data types and constraints used in the JHI in terms of root-mean squared error (RMSE). For JHI Run 1, head, flux and TDEM data are fitted with an RMSE of around 1 for most realizations. In JHI Run 4 coupling constraints become stronger and the RMSE for the flux and TDEM data start to increase. The head data do not clearly show this behavior. The RMSE for the petrophysical coupling constraint shows a decrease for JHI Runs 4 and 7, whereas the RMSE of the geometric coupling constraint increases. The latter demonstrates the dominance Hydrol. Earth Syst. Sci., 17, 4043-4060, 2013 www.hydrol-earth-syst-sci.net/17/4043/2013/ rror parameter estimates for the second configuration of JHI and SHI runs using 250 parameter realizations and larger e rated "true" parameters. Blue dashed lines indicate parameter errors for JHI Run 1-7, where the cyan lines indicate the er errors for the SHI. of the petrophysical coupling constraint due to the employed weighting strategy and the high parameter sensitivity of r 1 that is subjected to this constraint.
Secondly, a SHI is applied to evaluate the performance of the JHI. The cyan lines in Fig. 3 show the parameter errors for the SHI. These results show a large reduction in parameter error for K clay and D clay compared to JHI Run 1. For parameter K clay this reduction of error is similar to JHI Runs 6 and 7. For D clay the SHI performs better compared to JHI Runs 6 and 7, indicated by the number of JHI realizations with an error larger than 15 %. Compared to these runs the geophysical parameter errors are generally smaller for the SHI. The last column in Fig. 4 lists the data fit for the SHI. As the inverted TDEM models of JHI Run 1 are used in the SHI, the histogram for the TDEM data is identical to that of the TDEM data in JHI Run 1. Head and flux data are fitted less well compared to JHI Run 1. The fit for both coupling constraints indicate a relatively strong petrophysical constraint.
Finally, a second configuration of JHI and SHI is tested in which a larger standard deviation (e corr ) was used to generate less correlated parameter realizations for K clay , D clay , r 1 and t 1 ; 0.1 for the petrophysical constraint and 0.5 for the geometric constraint. Figure 5 shows a reduction in parameter error for K clay compared to JHI Run 1 from about 100 to 60 %. The SHI resulted in a similar reduction. The improvement in K clay , however, is much smaller compared to the results in Fig. 3. Geophysical parameters r 1 and r 2 in JHI Runs 6 and 7, show worse estimates compared to JHI Runs 6 and 7 in Fig. 3.
The average computational burden associated with the inversion for a single realization was 94 (61 + 33) model calls for the SHI compared to 306 (153 + 153) model calls for the JHI. As the estimation of geophysical and groundwater model parameters is conducted simultaneously, the number of iterations in which geophysical and groundwater model parameters are updated are the same, which is not the case in a SHI. This will result in a larger computational burden for the JHI.

Example 2: case study Risby landfill
As second example we consider a steady-state, real-world groundwater model for Risby landfill located in Denmark, to which we refer as the Risby model. This model was developed by Christensen and Balicki (2010) to characterize the hydrogeological interaction between a landfill, a local stream and a regional aquifer that is used for water supply. Christensen and Balicki (2010) provide a thorough description and discussion of the assumptions underlying the setup of this model and its results.
We investigate the application of a SHI and JHI to inform the groundwater model with ERT data that was collected near Risby landfill (Fig. 6). We first list the basic properties of the 1 Fig. 6 An aerial overview of Risby landfill, the ERT profile, parameter PP1 and available boreholes and hydrogeological observation 2 data at Risby landfill.
3 Fig. 6. An aerial overview of Risby landfill, the ERT profile, parameter PP 1 and available borehole and hydrogeological observation data at Risby landfill.
Risby area and the Risby groundwater model, after which we conduct a simple linear sensitivity analysis for the different hydrogeological parameters in the groundwater model, followed by the application of a SHI and JHI to inform the groundwater model with the ERT data.

Description of Risby landfill
An extensive historical overview of Risby landfill was provided by Thomsen et al. (2011). Figure 6 lists the key features of the study area, which are a landfill and a small brook called Nybølle stream. The geological setting of Risby landfill (Højberg et al., 2008;CarlBro, 1988) comprises pre-Quaternary limestone bedrock overlain by Quaternary glacial deposits. The pre-Quaternary limestone surface is located between −10 and +5 m a.m.s.l., corresponding to 20-30 m below the natural terrain surface. The Quaternary glacial deposits mainly consist of clay till, but intercalated sand lenses and sand layers are common. The sandy deposits range in thickness from a few centimeters to several meters. Figure 7a shows the horizontal grid discretization that is used to simulate groundwater levels near Risby landfill. The grid cell size employed in the groundwater model is 50 m by 50 m. Near the landfill a smaller cell size of 12.5 m by 12.5 m is employed. For the geological setup, 5 continuous layers were chosen, where the 4 upper layers represent the sand and clay layers of the glacial clay till and the lowest layer represents the field-scale limestone aquifer. The top layer of the model, with its bottom elevation fixed at +15 m a.m.s.l. was subdivided in 3 zones, which represent the extent of the upper sandy and clayey deposits together with the delineation of the northern part of the landfill (Fig. 7a). Boundary conditions applied in the Risby model are shown in Fig. 7b and consist of constant heads, derived from a regional groundwater model, referred to as the GEUSmodel (Højberg et al., 2008). The limestone was assumed to be impermeable at a level of −50 m a.m.s.l. and a no flow boundary was therefore assigned. The boundaries for the top layer and the remaining two clay layers were also set as no flow boundaries. The symbols Q GEUS , H GEUS and R GEUS indicate the specified flux, constant head values and recharge, which were extracted from the regional GEUSmodel. Boundaries for the limestone were set as constant head boundaries with a hydraulic head equal to 14.9 m. The isopotential used, was the average simulated head in the limestone for the period 2001-2005 (Højberg et al., 2008). Boundaries for the sand layer were prescribed flux boundaries. A flux of 7.2 × 10 −6 m 3 s −1 was applied for all cells along the boundary.

Groundwater model
In Christensen and Balicki (2010) the Risby model was calibrated using 6 parameters listed in Table 3, representing a uniform hydraulic conductivity for every geological layer, except for the uppermost layer, which consists of 3 separate zones and the bottom clay layer for which the hydraulic conductivity was fixed. The observation data comprised 34 head measurements and 4 flux measurements (Fig. 6).

ERT data
The landfill and its surroundings were mapped using various geoelectrical profiles for which ERT and induced polarization data (Slater, 2007) were collected in order to delineate the landfill, sand pockets and the thickness of the glacial deposits overlying the limestone aquifer (Gazoty et al., unpublished data). To demonstrate the SHI and JHI, we used the data associated with one of these ERT profiles north of the landfill; the location of the profile is shown in Fig. 6. Figure 8a shows the inverted resistivity model for the ERT profile using a few-layer, laterally constrained inversion (LCI) approach as discussed in Sect. 2.1. This ERT profile consists of 37 1-D resistivity models with 3 layers and is orientated west-east (model number 0 marks the western point). The parameters estimated for each of the 37 resistivity models (5 m spaced) comprise 3 layer resistivities (r 1 , r 2 and r 3 ) and 2 layer thicknesses (t 1 and t 2 ). Lateral constraints were used with a weight factor of 1.2 for the layer depths (C Rh ) and a weight factor of 1.2 for the resistivities between neighboring resistivity models. These weight factors are described in Auken and Christiansen (2004) and their value is subjectively determined and based on common practice ranges suggested in HGG (2008).
At the location of the ERT profile, boreholes showed a depression in the limestone surface down to ca. −10 m a.m.s.l. This depression has been interpreted as a buried paleovalley in the pre-Quaternary landscape and its shape is not well captured with the available boreholes. Another characteristic are relatively thick sand layers at the eastern part of Risby landfill.
In Fig. 8a  1 Fig. 8 Inverted ERT model obtained after a separate geophysical inversion (a) and using the J 2 parameter uncertainty analysis expressed by their standard deviation relative to the parameter 3 colored) and undetermined parameters (light colored) for the separate geophysical inversion ( 4 5 Fig. 8. Inverted ERT model obtained after a separate geophysical inversion (a) and using the JHI with e c = 0.2 (b) together with a parameter uncertainty analysis expressed by their standard deviation relative to the parameter estimate. A gray scale marks well (dark colored) and undetermined parameters (light colored) for the separate geophysical inversion (c) and a JHI with e c = 0.2 (d).
high electrical resistivities of about 50-80 m recorded at the eastern part of the profile (model numbers 15-37). The top layer with a resistivity of ca. 10 m is more pronounced at the western part of the profile (model numbers 1-10), indicating predominantly clayey deposits. The presence of the landfill and an associated leachate plume might slightly affect this estimated resistivity. Leachate migration is not considered in this study because the discretization of the groundwater model is insufficient to accurately simulate this process  (Milosevic et al., 2012). Figure 8c shows the uncertainty associated with the parameters that are estimated in the ERT model, expressed by their standard deviation as a percentage of the parameter estimate. This parameter uncertainty analysis included all the information provided by the data and parameter constraints. Note light colors in Fig. 8c indicate relatively poorly resolved parameters, e.g., r 1 , r 2 and t 1 for models 1-10.

Informing the Risby model with ERT data
As mentioned before, 6 parameters are estimated in the original Risby model (Christensen and Balicki, 2010), which are listed in Table 3. For these parameters a local, linear sensitivity analysis (Fig. 9) is conducted using PEST (Doherty, 2010). This analysis shows that the hydraulic conductivity pertaining to the clay layer (K clay ) is the most sensitive parameter.
To improve the estimate of K clay a petrophysical relationship is applied, which is used in Eqs. (25) and (26). An expected value of 9 is used, as clay till has an approximate hydraulic conductivity of 10 −8 m s −1 (Fredericia, 1990;CarlBro, 1988) and an electrical resistivity of about 10 1 m (Kirsch, 2006). This relationship implies a higher electrical resistivity is accompanied by a smaller clay content, which, in turn, results in a higher hydraulic conductivity. r 1 and r 2 in resistivity model numbers 1-10 are coupled to the estimation of K clay , as the area eastern part of the ERT profile (model numbers 15-37) contains large sandy deposits embedded in the clay. As we are only using a 3 layer resistivity model the average electrical resistivity in this part of the domain would not reflect the resistivity of the clay appropriately.
As the ERT model also informs about the depth to the limestone, we introduce an additional parameter (PP 1 ) in the groundwater model representing the top elevation of the limestone. PP 1 represents a single pilot point (Certes and Demarsily, 1991) used to interpolate the elevation of the limestone surface together with the available borehole information. The location of PP 1 , which is shown in Fig. 6, is picked as the depression of the limestone surface, occurring at the northeastern part of the landfill, is not well characterized. As expected, the calculated sensitivity, based on Hill (1998), of this parameter is very small with respect to the hydrogeological observations (Fig. 9). To demonstrate the effect of geometric coupling we use parameter PP 1 in the inversion process. Parameters t 1 and t 2 in model numbers 14, 15 and 16 are coupled to the estimation of PP 1 .

SHI
The SHI starts with the estimated geophysical model shown in Fig. 8a. The scale of the individual 1-D resistivity models comprised by the ERT model is rather small (electrode spacing of 5 m) compared to the grid cell size of 12.5 m used in the groundwater model. For this purpose we have chosen to constrain K clay with the average electrical resistivity estimates, r 1 and r 2 , pertaining to resistivity model numbers 1-10. To constrain the estimation of PP 1 we use the average sum of t 1 and t 2 pertaining to resistivity model numbers 14, 15 and 16. The weights associated with the constraints were based on the standard deviations of the geophysical parameter estimates calculated using Eq. (7).

JHI
We also apply a JHI for the Risby model to estimate r 1 , r 2 and K clay using the petrophysical relationship described in Sect. 4.2. For the estimation of the depth to the limestone we introduce a geometric coupling constraint between parameters PP 1 , t 1 and t 2 . The petrophysical coupling constraint is used for resistivity models 1-10, the geometric constraint for resistivity models 14, 15 and 16.

Results
The last column in Table 3 shows the parameter estimation results for a separate inversion of both the geophysical and the groundwater model. Most of the parameters in the groundwater model are estimated with relatively small posterior standard deviation. When performing a SHI (Table 3, column 2), the decrease in parameter uncertainty is negligible except for K lime and K risbyn . Parameter estimates remain similar to the separate inversion, which is likely caused by the large standard deviations associated with the geophysical parameters that are coupled with the groundwater model. In Fig. 8c these parameters show a relatively large standard deviation. As we used this standard deviation to determine the weight of the constraints in the SHI, the constraint might be too weak to affect the estimation of the groundwater model parameters significantly. Figure 10 shows the parameter estimates and 68 %confidence intervals for the JHI, when using different weight values for the coupling constraints (e c ). The parameter estimates for K clay and r 2 are affected when the weight of the petrophysical relationship is increased by setting the acceptable error e c to a smaller value. The geometric constraint between PP 1 , t 1 and t 2 does not have a big impact on the Hydrol. Earth Syst. Sci., 17, 4043-4060, 2013 www.hydrol-earth-syst-sci.net/17/4043/2013/ 1 Fig. 10 Parameter estimates (black straight line) and confidence bounds (red dashed lines) for different values of e c when performing 2 JHI using a petrophysical relationship between K clay , r 1 and r 2 and a geometrical constraint between parameters PP 1 and t 1 and t 2 . The 3 confidence bounds represent the parameter estimate ± 2 standard deviations. 4 5 6 Fig. 10. Parameter estimates (black straight line) and confidence bounds (red dashed lines) for different values of e c when performing a JHI using a petrophysical relationship between K clay , r 1 and r 2 and a geometrical constraint between parameters PP 1 and t 1 and t 2 . The confidence bounds represent the parameter estimate ± 2 standard deviations. estimated values of the geophysical parameters. However the estimate of PP 1 does approximate the geophysical model better when the constraint is given more weight. The average depth to the limestone in the ERT model is about 25 m (t 1 + t 2 ). In the groundwater model, this depth is estimated to be 28.3 m ± 0.5 and 28.0 m ± 1.0 m using a separate inversion and a SHI, respectively. In the JHI this estimate becomes ca. 26.6 m ± 0.6 m. Table 3 shows that standard deviations of the groundwater model parameters for the JHI are almost equivalent compared to the SHI, but smaller compared with a separate inversion.
The main advantage of the JHI is seen from the estimated values for the geophysical parameters that are allowed to change in the JHI. Geophysical layer thicknesses, t 1 and t 2 , decrease slightly compared with the SHI, while electrical resistivity r 2 shows a more significant change. Figure 8b is the inverted ERT model using the JHI with an e c of 0.2. Compared with the geophysical inversion results in Fig. 8a the estimated resistivity of layer 2 dropped from an average of 75 to ca. 30 m for resistivity models 1-10. These are the models for which electrical resistivities r 1 and r 2 were coupled to K clay in the groundwater model. Figure 8d shows the standard deviations associated with the estimated geophysical model obtained with the JHI. The standard deviation of parameter r 2 indicates it is not well determined using the JHI as was the case in the separate geophysical inversion. r 1 is determined with an approximate standard deviation of 10 %. However, Fig. 8d shows t 1 is less well resolved for those model numbers where the petrophysical relationship is applied. The geometric coupling constraint does not show any effect on the estimated geophysical models in Fig. 8. Table 3 lists the RMSE with respect to the geophysical and hydrogeological observations (respectively ϕ g and ϕ h ), which was smaller than 1 for all simulations. No significant increase in data fit was noted, except a slightly higher ϕ h for the JHI. Increasing the weight of the coupling constraints (by decreasing e c ) or increasing the number of coupling constraints, will ultimately result in an increase in ϕ g and ϕ h , as the geophysical and groundwater data will pull parameters in different directions.
The last entry in Table 3 is the amount of model runs needed to perform the different inversion types. The JHI required about twice as many geophysical and groundwater model runs compared to the separate inversion and ca. 3 times as many groundwater model runs compared with the SHI.
This study tested a SHI and a new type of JHI for a groundwater model and different types of geophysical data. The JHI estimated geophysical and groundwater parameters simultaneously, employing coupling constraints acting as additional regularization terms to exploit potential correlation between geophysical and hydrogeological properties that can be based on established petrophysical relationships. The SHI employed similar coupling constraints, but included an independent geophysical inversion. The weight of the SHI coupling constraints was based on geophysical parameter resolution.
Both the SHI and JHI approaches can provide consistent inversion frameworks and offer a high level of flexibility when coupling groundwater and geophysical models because 1. only selected geophysical model parameters can be coupled to groundwater model parameters, 2. confidence associated with the hydrological interpretation of a geophysical model can be tuned using different weights for the employed coupling constraints, 3. scale issues can be overcome by coupling several geophysical parameters to hydrological parameters or vice versa, 4. SHI and JHI can be applied for various combinations of geophysical methods and groundwater models, and 5. SHI and JHI can be used with other types of optimization methods (e.g., Markov-Chain Monte Carlo methods) by adding an additional coupling constraint component to the objective function that is minimized.
Furthermore, the JHI and SHI are consistent with state-ofthe-art inversion techniques used for groundwater models, resistivity and airborne electromagnetic data. For a synthetic study, comprising a cross-sectional groundwater model and TDEM data, a JHI and SHI resulted in improved parameter estimates and a reduction in parameter uncertainty in comparison with a groundwater model that is not informed with TDEM data. Groundwater parameter estimates using a JHI did not improve compared with a SHI and resulted in slightly worse parameter estimates for the geophysical model when using large weights for the coupling constraints. A second configuration of the synthetic study, incorporating lower quality (petro)physical relationships between geophysical and groundwater parameters resulted in decreasing performances for both the SHI and JHI. The SHI performed slightly better compared to the JHI based on the geophysical parameter estimates and geophysical data misfit. In contrast to the JHI, the SHI overestimated the level of correlation between geophysical and groundwater parameters. To avoid overestimating model coupling strength in a SHI (which can result in an underestimation of parameter uncertainty), weighting strategies for parameter coupling constraints should be based on that element (parameter resolution or petrophysical relationship) that incorporates the largest error.
For the case of a real-world, field-scale groundwater model and an ERT section, parameter uncertainty was significantly decreased for two parameters in the groundwater model using both a JHI and SHI. The JHI resulted in different parameter estimates for both the groundwater and the geophysical model, honoring the imposed coupling constraints. Parameter uncertainty was not reduced in comparison with a SHI.
For the cases investigated in this paper the SHI proves to be more useful based on analyses of parameter estimates and data fit. In addition, the JHI requires a 2-3 times larger computational burden and is relatively difficult to implement. The JHI might still be useful when groundwater and geophysical models can mutually benefit from differences in parameter resolution. For coupling geophysical models with field-scale or regional groundwater models, such situation is not likely to occur as the groundwater models are relatively more prone to conceptual errors and limited observation data. Finally, when planning hydrogeophysical surveys and modeling, parameter sensitivity studies are of crucial importance to explore parameters that need to be determined, given targeted groundwater model predictions, and to determine whether parameter resolution in geophysical models provides opportunities to constrain these parameters.