In many subsurface industrial applications, fluids are injected into or withdrawn from a geologic formation. It is of practical interest to quantify precisely where, when, and by how much the injected fluid alters the state of the subsurface. Routine geophysical monitoring of such processes attempts to image the way that geophysical properties, such as seismic velocities or electrical conductivity, change through time and space and to then make qualitative inferences as to where the injected fluid has migrated. The more rigorous formulation of the time-lapse geophysical inverse problem forecasts how the subsurface evolves during the course of a fluid-injection application. Using time-lapse geophysical signals as the data to be matched, the model unknowns to be estimated are the multiphysics forward-modeling parameters controlling the fluid-injection process. Properly reproducing the geophysical signature of the flow process, subsequent simulations can predict the fluid migration and alteration in the subsurface. The dynamic nature of fluid-injection processes renders imaging problems more complex than conventional geophysical imaging for static targets. This work intents to clarify the related hydrogeophysical parameter estimation concepts.

There are many scenarios in which the earth’s subsurface is being altered by anthropogenic activity. Our focus will be on processes that involve the injection of fluid into the subsurface such as that which occurs during wastewater storage, hydraulic-fracture generation, CO2 sequestration, enhanced geothermal-energy generation, enhanced oil recovery, groundwater remediation, and underground gas or liquid storage. Any of these applications involve (hydrologic) state changes in some subsurface system due to fluid injection. Understanding and predicting its impact requires quantitative estimates of not only where the fluids went once they were injected but also how they altered the state of the subsurface. As will be outlined in detail, these concepts and overall goals are distinct from those underlying time-lapse geophysical monitoring approaches such as crosswell seismic (e.g., Landrø and Stammeijer, 2004; Daley et al., 2008; Marchesini et al., 2017), seismic coda monitoring (e.g., Kanu et al., 2014; Obermann et al., 2016), crosswell electromagnetics (EM) (e.g., Binley et al., 2001; Day-Lewis et al., 2003; Marsala et al., 2008), and crosswell electrical resistivity tomography (ERT) (e.g., Daily et al., 1992; Bergmann et al., 2012).

In this tutorial, we elaborate on the concepts around estimation of parameters that control the migration of injected fluids. Despite our specific application, we suggest a broader point of view, called “imaging a subsurface process.” The main distinction of the process-imaging approach presented here is in the definition of the forward model that simulates how time-lapse geophysical signals vary during the course of fluid injection. In conventional formulations of geophysical-monitoring inverse problems, the unknowns are the time-varying geophysical properties such as seismic velocities, electrical conductivity, or the dielectric constant, in which the forward model only simulates the corresponding geophysical signals (seismic traces, electrical voltages, etc.). The process-imaging formulation entails modeling of geophysical signals as one component of a larger coupled multiphysics forward model.

Our reasoning for the process-based perspective is that geophysical observations have found their way into a range of subdisciplines that involve various flow processes, an early one being reservoir engineering (e.g., Nolen-Hoeksema, 1990; Pagano et al., 2000), followed by the vast field of hydrogeophysics (e.g., Rubin and Hubbard, 2005a; Vereecken et al., 2006; Binley et al., 2015); soil moisture studies (e.g., Samouëlian et al., 2005), biogeophysics (e.g., Atekwana and Slater, 2009), agricultural geophysics (e.g., Allred et al., 2010), and geothermal exploration (e.g., Bromley, 2018).

The trend of geophysical inquiry of hydrologic systems has given rise to a variety of deterministic and stochastic hydrogeophysical parameter estimation schemes, along with a variety of classification attempts (Linde et al., 2006a; Ferré et al., 2009; Hinnell et al., 2010; Herckenrath et al., 2013; Liang et al., 2014; Camporese et al., 2015; Linde and Doetsch, 2016). Here, we will follow the simple distinction between (fully) coupled and uncoupled hydrogeophysical inversion methods (e.g., Hinnell et al., 2010) because we will show joint inversion examples that are closely related to the coupled method. Further explanation and more related methodological reviews are given below.

Our didactic treatment of the subject also tries to bridge a possible conceptual gap that may exist when practitioners transition from conventional (often large-scale) geophysical inversions for static targets — common in exploration and reconnaissance problems — to imaging of dynamic flow systems. Hence, we will first lay out the overall concepts in a way that is applicable to any imagined fluid-injection process as constrained by geophysical data. Afterwards, we make the ideas definite by giving a particular example of a multiphysics forward model, namely, the injection of salt water into a permeable formation with time-lapse electrical data gathers used to image the process. The last portion of the paper presents a numerical approach for performing multiphysics inverse modeling for process-driving parameters. Again, mainly for didactic purposes, these demonstrations attempt to be distinct from the hydrologic literature by using techniques that are more common to conventional deterministic geophysical imaging on rectangular parameter grids.

Finally, in discussing this particular forward-modeling application, we attempt to be self-contained. Excursions into related hydrogeophysical aspects are interspersed throughout on a contextual basis, assuming some textbook-level background of basic hydrogeology (e.g., Helmig, 1997; Hiscock and Bense, 2014). In discussing the inverse problem, we also assume that the reader has some familiarity with least-squares inverse modeling concepts as presented, for example, by Menke (1984) and Tarantola (2005). Vasco and Datta-Gupta (2016) cover imaging and subsurface fluid flow principles with a focus on trajectory-based flow modeling.

By “imaging the process” of fluid injection, we mean that we have the ability to simulate numerically how the injected fluid, along with the solute and heat it carries, enters the subsurface, migrates from the borehole, and changes pertinent subsurface properties. To be successful, we must have confidence that such simulations represent, to a reasonable approximation, what actually has occurred in the subsurface. The ability to simulate a process and have the simulation be realistic has two equally important components.

First, we need forward models of the fluid-injection process that appropriately represent the pertinent physics and chemistry that take place during the fluid-induced subsurface alterations. They typically come in the form of sets of coupled partial differential equations (PDEs) that are solved numerically. Such multiphysics forward modeling is also referred to as THMC (for thermal, hydraulic, mechanical, and chemical) models (e.g., Taron et al., 2009). If necessary, models are also provided for how the coefficients of the PDEs — or material properties such as permeability and porosity — evolve over time. Bromley (2018) alludes to this aspect while discussing the role of geophysical monitoring of THMC processes. The dependent variables of the PDEs will be called here the “primary fields”; they are fields such as fluid pressure, Darcy velocity, chemical concentration, temperature, stress, particle displacement (or particle velocity), and electric fields.

The forward models involve voxel resolution that is appropriate to the fluid-injection application and always use voxels that are much larger than the grain sizes of the material. Hence, simulations of the fluid-injection process use a macroscopic “porous-continuum” description of the physics and chemistry. Much of the basic research surrounding the overall topic of imaging subsurface processes is in improving the porous-continuum forward models that describe the fluid-injection process, including the associated subsurface property changes.

Second, in order for such forward modeling to be a useful representation of reality, we must know the initial conditions of all fields in the forward model throughout the portion of the subsurface that will be affected by the fluid injection. These initial conditions include not only the spatial distribution of the primary fields of the process being simulated but also the spatial distribution of the various material properties that are the coefficients in the PDEs governing the fluid-injection process. These material properties are quantities such as Darcy permeability, porosity, thermal conductivity, thermal capacity, solute dispersivity, and fluid-storage capacity.

We will refer to the material properties pertinent to the forward modeling of the fluid injection process as hydrologic properties. In contrast, geophysical (material) properties refer to the forward-modeling PDE coefficients that define how time-lapse geophysical signals vary during the course of fluid injection. Examples are electrical conductivity or elastic moduli. In addition, to distinguish from geophysical forward-modeling processes, we will consistently use the adjective “hydrologic” when referring to any aspect related to the fluid injection process.

The process-imaging inverse problem is to obtain all relevant initial primary fields and material property fields, as well as any other process-relevant parameters that are required to simulate how physical properties evolve through time. All three groups define the set of input parameters that need to be calibrated through some parameter estimation scheme, i.e., inversion procedure.

Once all input parameters have been properly calibrated through the inverse analysis, the forward model’s output matches given observations. We can then say that we are imaging the subsurface process given the observational database. Note that in addition to time-lapse geophysical signals, the observations can include measurements of the hydrologic primary fields (e.g., fluid pressure, temperature, or chemical concentration measurements in boreholes), as well as samples of material properties from well logs or core data.

Formulating the process-imaging inverse problem

Geophysical imaging commonly refers to the procedure of analyzing recorded geophysical signals to obtain maps of the spatial distribution of geophysical properties such as seismic reflectivity, seismic velocity, or electrical conductivity. If time-lapse collection of geophysical signals is taking place, one can monitor how the geophysical properties vary through time in the subsurface. One then qualitatively associates a changing geophysical property at a voxel with the arrival of injected fluid or altered stress at that voxel. In the seismic industry, such time-lapse 3D imaging over petroleum reservoirs is also known as 4D monitoring (e.g., Fanchi et al., 1999).

The coupled (hydrologic and geophysical) forward problem simulates the fluid-injection process in conjunction with induced time-lapse changes of seismic or electrical properties and corresponding changes of artificially sourced geophysical signals. Geophysical data gathers serve to constrain the unknowns in the corresponding parameter estimation process. Therefore, we have one large coupled multiphysics forward model that is predicting, among many things, how measurable geophysical signals are changing as the fluid invades the subsurface.

The overall inverse workflow can then be summarized as (we assume the geophysical data to be the only given observations):

  1. 1)

    Forward-model the entire fluid-injection process from the beginning of fluid injection through the latest time-lapse geophysical data. The forward modeling includes how the geophysical properties and geophysical time-lapse signals are being altered by the fluid injection. Rock-physics models for how the pertinent material properties are being altered by the fluid injection are part of this multiphysics forward model.

  2. 2)

    Evaluate the difference between the recorded and simulated time-lapse geophysical signals. If a satisfactory data fit is achieved, stop. Otherwise, adjust the forward-modeling input parameters with an influence on the observations.

  3. 3)

    Return to step 1 and repeat until an optimal estimate of the input parameters is obtained throughout the region being affected by the fluid injection at the current time-lapse geophysical data gather.

This scheme follows the general framework of the coupled hydrogeophysical inversion approach described by numerous authors, which we want to recount in the following section.

Excursion: Coupled versus uncoupled hydrogeophysical inversion

Step 1 in the parameter estimation workflow outlined above tightly couples the geophysical signal simulation to the fluid-injection forward model. One could say that the geophysical (material) properties from which data are predicted are constantly being updated internally through the fluid injection process. This workflow is also known as (fully) coupled hydrogeophysical inversion, with its counterpart often referred to as an uncoupled (also decoupled or sequential) approach. In uncoupled approaches, the geophysical properties pertinent to a fluid-injection process are obtained via stand-alone geophysical inversion.

Numerous uncoupled approaches have been reported in the literature (e.g., Hyndman and Gorelick, 1996; Binley et al., 2002; Kemna et al., 2002; Lambot et al., 2004; Singha and Gorelick, 2005; Vanderborght et al., 2005; Linde et al., 2006b). Their commonality is that standalone geophysical data inversions, whether with one or more data types, are part of the workflow. Resulting tomograms can be converted to hydrologic properties and/or spatial plume delineations. The linkage is typically given by established rock-physics (or petrophysical, where we use both terms interchangeably) relationships, or through structural similarity constraints. Subsequent hydrologic model calibrations use the converted tomograms in the form of a priori information and/or hydrologic proxy data.

Few studies compare the uncoupled against the coupled approach (Sicilia and Moysey, 2007; Hinnell et al., 2010; Camporese et al., 2015). The main drawback of coupled approaches is a bias from incorrect rock-physics assumptions, whereas the decoupling allows for alternatives through more forgiving ways such as structural joint inversion (e.g., Linde et al., 2006a; Lochbühler et al., 2013). The main argument in favor of coupled approaches is that they enforce a quasi regularization through the conservation and transport laws of the fluid-injection process itself, thus avoiding the need for other regularization such as smoothing (Hinnell et al., 2010; Camporese et al., 2015). In other words, the uncoupled geophysical inversion is not regularized by the physics posed by the fluid-injection process, thus needing regularization to mitigate its ill-posed nature. It is typically given through smoothing or other a priori information, which may not always be based on physical principles or site-specific characteristics and thus may introduce unwanted biases. Therefore, diffuse tomograms due to excessive smoothing, or other inversion artifacts, can adversely affect the final hydrologic image of uncoupled inversions (e.g., Hinnell et al., 2010).

Our didactic purpose calls for some clarification on the regularization aspect because the argument that coupled inversion overcomes the ill-posed nature of smoothness-constrained geophysical inverse problems may lead to confusion. The reason is that this argument is somewhat inessential because the inversion for geophysical material properties is absent in a coupled inversion for the underlying (process-defining) hydrologic material properties. In the presence of insufficient data, every overparameterized inverse problem becomes ill-posed. Whether imaging a spatially complex hydrologic state (e.g., tracer concentration) through uncoupled resistivity inversion or imaging a spatially complex distribution of underlying hydrologic properties that control the state (such as permeability) through coupled inversion, both cases will require regularization. In other words, although the evolution of resistivity is controlled by the physics of the flow system, the distribution of permeability is not. A dedicated numerical example below will support this clarification.

Thus far, the high computational demands have rendered the regularization in coupled inversions less of a pressing need because most demonstrations have involved only few parameters. They describe homogeneous models (Rucker and Ferré, 2004; Sicilia and Moysey, 2007; Lehikoinena et al., 2009), horizontally layered (1D) models (Looms et al., 2008; Hinnell et al., 2010; Irving and Singha, 2010; Kowalsky et al., 2011; Scholer et al., 2011; Busch et al., 2013; Tran et al., 2014; Vilhelmsen et al., 2014), or preset facies-like or dike structures (Huisman et al., 2010; Rings et al., 2010). Geostatistical methods permit more spatial variability through structural and stratigraphic information (Kowalsky et al., 2004, 2005; Finsterle and Kowalsky, 2008; Jardani et al., 2013; Ahmed et al., 2014, 2016; Camporese et al., 2015), also using a handful of forward-modeling input parameters numbering several orders below that of typical large-scale geophysical 3D imaging problems.

Calculating large parameter sensitivity matrices for Gauss-Newton style coupled inversions becomes computationally feasible when certain simplifying assumptions can help accelerate the forward modeling for either the fluid-injection process or the geophysical signal evolution process. Faster flow simulation can, for example, be achieved through trajectory-based modeling, sometimes also referred to as streamline simulation, which basically resembles ray-based seismic imaging (Vasco et al., 2004). On the other side, ERT sensitivity calculations can be somewhat condensed through the use of temporal moments of electrical potential perturbations, which essentially lets one invert mean arrival times of ERT signals, instead of full time series of voltages (Pollock and Cirpka, 2010; Pollock and Cirpka, 2012). We note that the latter three references present coupled inversion studies in which the objective functions have regularizing terms that penalize large spatial parameter fluctuations.

Fluid-injection effects on material properties

The above inversion workflow is based on the notion that the injected fluid is altering the state of the subsurface, which in turn alters a geologic media’s material properties and thus geophysical signals that are being recorded during the course of the application. The fluids being injected into the subsurface are out of equilibrium with the in situ fluids. They will necessarily be at a higher fluid pressure, may be at a different temperature, and will typically be out of chemical equilibrium, that is, the solutes and water will have different chemical potentials compared to the in situ fluids.

There are four main ways that injected fluids alter the material properties of rock. First, fluid substitution is the process of an in situ fluid being replaced by the injected fluid. The latter usually has different physical properties, caused by differences in chemistry, temperature, and pressure. Focusing specifically on seismic rock properties, the tutorial by Smith et al. (2003) details how to quantify fluid substitution effects. The second way involves stress-induced changes to rock properties. Here, the injected fluid provokes deformation and stress change, sometimes at a considerable distance from where fluid substitution is occurring. As a consequence, preexisting cracks and fractures can close or open, further altering the mechanical and transport properties (e.g., Pride et al., 2017b). Third, fluid injection may induce rock damage. Corresponding stress changes can provoke the creation of new grain-scale cracks, larger scale fractures, or slip on preexisting faults, which can further alter the elastic and transport properties. Finally, chemical nonequilibrium introduced by an injected fluid may infer precipitation or dissolution of minerals, also altering the elastic and transport properties (e.g., Kang et al., 2003).

The following specific example involves a higher salinity solution that is injected from a well into a heterogeneous geologic formation initially saturated with a lower salinity solution. We assume that a background flow field exists in addition to the flow created by the injection of the saline solution. To keep the example simple, we assume that the porous continuum is isotropic, though heterogeneous, and that there is a single in situ fluid that fully saturates the pore space. The time-lapse geophysical data acquisition used in this example is an ERT survey carried out across a group of monitoring wells.

Modeling the fluid injection and associated salt dispersion

In the following particular example, fluid substitution is the sole mechanism for provoking changes to the pertinent geophysical property, which is electrical conductivity.

The saline solution, injected and in situ, is assumed to be an NaCl solution in which the Na+ and Cl ions are completely dissociated in solution, a so-called strong electrolyte. We quantify the salt concentration in each voxel of porous material using the mass ratio c, which is defined as the total mass of Na+ and Cl ions in the voxel divided by the total mass of solution (Na+, Cl, and H2O) in the voxel.

The solute mass balance in a porous material is given by (e.g., a recent first principle derivation is given by Pride et al., 2017a)
where q is the Darcy filtration velocity, ρf is the total solution mass density of the fluid, D is a dispersion tensor, and ϕ is the porosity so that q/ϕ is the average speed that the fluid is moving through the pores relative to the solid grains. The advective derivative on the left side largely represents solute accumulation due to the average advection of solute, that is,
with ·q corresponding to the total fluid density changes and ·(cq) corresponding to the solute accumulation due to the average advection. The first term on the right side of equation 1 represents solute accumulation due to the diffusive and dispersive solute transport. The source (second) term on the right side corresponds to solution of concentration cs being injected at the volumetric rate Qs (m3/s) at the point rs in the subsurface. Note that this source term can have a nonconstant time history.
For isotropic materials, the second-order dispersion tensor D of equation 1 is often modeled as (e.g., Scheidegger, 1961)
where x^q is a unit vector in the direction of fluid flow defined as x^q=q/|q| and F is the formation factor, which is commonly modeled using Archie’s (1942) law as F=ϕm, where the so-called cementation exponent m depends on the grain shape. In what follows, however, we elect not to use Archie’s law and will treat the formation factor as an independent unknown property to be inverted for as it influences the electrical conductivity of the rocks.

In equation 3, γl and γt are called the longitudinal and transverse dispersivities that have units of length and have been measured in a host of laboratory experiments as being proportional to grain diameter d (e.g., Delgado, 2007), although they are often thought to increase with the increasing size of the plume (e.g., Gelhar et al., 1992).

Furthermore, in equation 3, the solute molecular diffusivity Dm is given by (e.g., Cussler, 2009)
where z± are the ionic valences that are the number and sign of fundamental charges on each ion (so z+=1 and z=1 for NaCl), ν± are the numbers of cations and anions in each salt molecule (so ν+=ν=1 for NaCl), and kBTb± are the ionic diffusivities of the individual Na+ and Cl ions. The form of equation 4 is derived from the requirement of charge neutrality during diffusion. Its other quantities are kB, Boltzmann’s constant; T, the temperature in Kelvin; and b±, the mobilities of the cations and anions as given by the Einstein-Stokes relation:
Here, the effective ion radius for Na+ is R+=1.63×1010 m and Cl is R=1.07×1010 m, while ηf is the solution viscosity.

To address the importance of diffusion to the evolution of the salinity plume, Appendix  A contains further details on a means of quantifying the ratio of advection to diffusion through a local Péclet number.

For the flow in the porous material, we combine the fluid mass balance with the compressibility law for the fluid in the porous material to obtain the rule for fluid-pressure change (see Pride, 2005)
which is complemented by Darcy’s law for the fluid flux through the porous material,
It is the heterogeneity of the permeability field across the modeling voxels that dominates the shape and evolution of the solute plume c(r,t) of interest. In the above, P is the fluid pressure, k is the permeability, g is the acceleration of gravity, Kf is the fluid bulk modulus, and Qs(t) is the volumetric rate that electrolyte is being injected into the subsurface at point rs. These laws for fluid flow assume that the fluid pressure changes are sufficiently small. Hence, permeability and porosity are not changing; that is, poroelastic deformation can be neglected. Otherwise, if the fluid injection is responsible for significant poroelastic deformation, then the laws of poroelasticity must also be included in the forward model.

The above transport equations are solved numerically with a finite-volume approach realized in the flow and transport simulator TOUGH2 (Pruess, 2004). The fluid properties that are influencing the above transport are the fluid density ρf, the fluid viscosity ηf, and the fluid bulk modulus Kf. The fluid properties ρf, ηf, and Kf vary with concentration c, fluid pressure P, and temperature T. TOUGH2 includes empirical relations for such state dependence of the fluid properties. In the present modeling example, the fluid injected into the subsurface is taken to be the same temperature as the subsurface, which we assume has the uniform value of T = 25°C.

Modeling the geophysical signals

The set of equations 17 describe the transport of the injected electrolyte through porous rocks, which has the effect of changing the rock’s electrical conductivity. To establish the link to geophysical measurements of this effect, we introduce the expression for the electrical conductivity of the electrolyte, which is given in SI units as
where μ± are the atomic masses of cations and anions (given in grams for one mole of these ions) so that ν+μ++νμ=58.44  gmol1 for NaCl, NA is Avogadro’s number (6.022×1023 objects mol1), ρf is the solution mass density, and e=1.60×1019  C is the fundamental charge. The factor of 1000 converts grams to kilograms so that σf is given in SI units. The ratio on the right side of equation 8 converts concentration c as expressed as a mass ratio to concentration as expressed in terms of numbers of NaCl molecules that went into unit volume of solution.
The time-lapse increase in electrolyte concentration c in a given porous voxel creates an increase in σf, further causing the bulk rock conductivity σR to vary through time according to the model
Equation 9 provides the geophysical material property input σR to our geophysical forward simulator.

The formation factor F in equation 9 depends only on the pore-space topology, and it is independent of the changing salinity field. If the framework of grains is highly deformable, the formation factor could vary with the pore pressure associated with the fluid-injection process. However, this effect is small for normal injection scenarios at depth (Pride et al., 2017b) and will be ignored here. Equation 9 is valid whenever electrical conduction along the surface of the grains of the rock is negligible. In this case, bulk ionic migration through the pores of the rock is the sole contributor to electrical conduction.

The geophysical signals due to changes in σR are recorded by means of a predefined array of ERT source and receiver electrodes, here given by a crosswell survey configuration. The source electrodes inject a current I (measured in Cs−1) into the geologic media, usually in a time-harmonic manner so that overvoltages at the sensor electrodes do not have sufficient time to build up. Such electrode impedances are associated with the conversion from ionic conduction in the earth to electronic conduction in the electrodes. However, the frequency of the applied current is sufficiently small that the associated EM diffusive skin depth is much greater than the depth of investigation so that the response is quasistatic in nature. The electric field in this case is expressed in terms of an electric potential as E=φ where the electric potential φ (the voltage) is the solution of the quasielectrostatic problem
with r+ being the electrode position where the possibly time-varying current I(t) is being injected and r being the position where the same amount of current is being extracted. Further details on the simulation of ERT data acquisition are given, for example, by Singha and Gorelick (2005), Vanderborght et al. (2005), and Pollock and Cirpka (2010).

Equation 10 is the final component of our multiphysics forward model in this example. One may note some conceptual similarity between the two interacting physical systems, one describing fluid injection (equation 1), while here (equation 10) it is an injection of electrical current. Both of their responses depend on their particular material properties.

Every time that we gather electrical data during the course of the fluid-injection experiment, we model the evolving rock conductivities due to fluid injection according to the above salinity model and its petrophysical link (equations 8 and 9). We then numerically solve for the voltages φ, which become data predictions to be compared against actual voltage measurements during the process-imaging workflow.

The fundamental difference to more traditional geophysical inversion, where one would solely estimate σR, is that the parameters adjusted in each subsurface voxel to match the voltage measurements are those that influence the transport of salinity, namely, the permeability k(r) and the formation factor F(r). Note that the formation factor is influencing the evolution of salinity through equation 3 and the rock’s electrical conductivity through equation 9.

The spatial variation of porosity ϕ(r) will always be much smaller in a percentage sense than the spatial variation in permeability k(r) and formation factor F(r). The places where porosity explicitly enters the above formulation, as a factor multiplying the time derivatives, are shown in equations 1 and 6. Note that the occurrences in equation 3 are neglected (Appendix  A). So even if the porosity goes through a small but sharp contrast at an interface, the porosity is not explicitly acted upon by a spatial derivative. As such, we will take the porosity to be a uniform constant in the inverse modeling, which is an approximation that can be tested.

Therefore, the spatial distributions of the material properties k(r) and F(r) represent the parameter space in the process-imaging inverse problem discussed next. All of the other fields are simulated through time and space according to the above multiphysics forward model. The initial condition for the fluid pressure prior to the start of electrolyte injection should also, formally, be a target of the inversion, but in the present example we assume that it is estimated prior to the tracer injection from well measurements. Similarly, the initial salinity field is assumed to be uniform throughout the subsurface and known prior to the start of fluid injection from well measurements. Due to such a priori information, the initial conditions of primary fields in the present model are assumed known.

Excursion: About petrophysical relations

Equation 9 constitutes the crucial coupling between two different physical systems, one describing the flow of injected fluid, and the other describing the flow of electric current. It is an essential element linking two different physical systems through a petrophysical relation between the evolving hydrologic primary fields and the evolving geophysical properties.

Most hydrogeophysical inverse modeling of ERT data assumes a linear relation between (electrical) rock conductivity and fluid conductivity (e.g., Singha and Gorelick, 2006a, 2006b), as is also assumed here. However, if conduction along the grain surface is allowed for in rocks with significant clay content and therefore a large grain surface area, the relation between σR and c (equations 8 and 9) becomes nonlinear. This can be addressed by additional terms in the petrophysical relation, which introduce additional material properties beyond the formation factor (e.g., Johnson et al., 1986; Pride, 1994).

For uncoupled inversions, alternatives have been suggested that replace potentially erroneous petrophysical assumptions by mere correlations. One approach replaces the petrophysical relationship by a correlation operator, the latter based on the assumption of a strong correlation between geophysical and hydrologic property changes (Johnson et al., 2009). Lochbühler et al. (2013) use spatial correlations between hydrologic and geophysical properties by enforcing structural similarity constraints.

To account for spatially changing petrophysical relationships, we choose to include petrophysical function parameters in the inversion, a common practice for coupled inversions (Hyndman et al., 2000; Kowalsky et al., 2005; Linde et al., 2006b, 2006c; Hinnell et al., 2010; Vilhelmsen et al., 2014).

The following synthetic data inversion examples involve the injection of a nonreactive saline tracer into a permeable formation that is bounded above and below by impermeable aquitards. Injection occurs from a single borehole over a period of 161 days with injection perforations located between a depth of 3.2 and 5.4 m below the upper aquitard.

Our goal is to simulate how the salinity within the heterogeneous confined aquifer evolves in space and time without having a priori knowledge of the heterogeneity of k(r) and F(r) present at the start of injection. A 2D model is considered here for didactic and computational simplicity. In addition, note that many tracer experiments of this kind involve only one well pair, making 2D models a first choice (e.g., Kemna et al., 2002; Vanderborght et al., 2005). The first inverse-modeling example considers a model with heterogeneous permeability and uniform formation factor. Afterward, we consider a second model in which permeability and formation factor are heterogeneous. The demonstration uses the inverse flow-modeling tool MPiTOUGH2 (Commer et al., 2014), which combines the inversion tool iTOUGH2 (Finsterle et al., 2017) with geophysical forward-modeling modules.

The geophysical model and data

The geophysical data portion is given by an ERT cross-borehole configuration. Figure 1 illustrates the flow domain that has dimensions of 3 m in height × 14 m in width. The in situ fluid is assumed to have an initial NaCl concentration of c=0.01, which corresponds to a molarity of approximately 0.1 M. Within the flow domain, the background electrical conductivity is altered over time due to the brine infiltration. Figure 1 shows a snapshot of the rock’s electrical conductivity due to elevated brine levels 19 days after injection began.

The geophysical data are given by three sets of ERT data, taken at 30, 70, and 110 days. The borehole array has 30 electrodes distributed over four wells. Each ERT data set contains 270 dipole-dipole configurations that are switched across the wells, generating 270 borehole-to-borehole ERT data points and totaling 810 points over the three data gathers. The electric potentials are modeled as having zero voltage along the domain boundaries.

The hydrologic model and data

The vertical plane of the flow domain is aligned with the direction of the background groundwater flow, and it spans 14 m in the horizontal direction and 3 m in the vertical direction (Figure 1). A total of 5098 grid blocks form the Cartesian forward-modeling domain. Grid intervals average 0.1 m in the zone between the injection and monitoring wells and 0.25 m near the left and right model boundaries. For the sake of this illustrative example, no fluid flow is assumed to occur across the top and bottom boundaries. Fixed pressure profiles are given to the vertical boundaries to impose a groundwater gradient of 56Pa/m in the lateral direction. Amended groundwater having c=0.1 (or a NaCl molarity of 1 M), which acts as an electrically conductive tracer, is injected uniformly throughout a well at x=3  m, indicated by the white line in Figure 1. Its density is higher than that of the resident groundwater due to the enhanced NaCl molarity of the injected electrolyte. Subsequent density-dependent effects are accounted for in the hydrologic simulations.

Figure 2 shows the actual heterogeneous permeability distribution, with values ranging from 4.77 E−13 to 1.90 E−10 m2. The heterogeneity is obtained through unconditional sequential Gaussian simulation (Deutsch and Journel, 1992) with an exponential semivariogram model for the logarithmic (log10) permeability with a sill of 0.26 and integral scales of 2.2 and 0.44 m in the horizontal and vertical directions, respectively. The permeability distribution influences how the injected brine spreads spatially over time.

In the following, we also refer to the field of brine concentration (the salinity field) as the C field. To see the effect of additional hydrologic data on the inversion, we assume that salinity measurements are made at four wells in the flow domain (Figure 2). These measurements are referred to as C data and are assumed to correspond to the spatially averaged salinity that is predicted by the flow forward model in the heterogeneous aquifer traversed by each monitoring well. Figure 2 shows the C field after 19 days of brine injection. Concentration data (C) are given as a time series of tracer (brine) concentration measured over flow time. Our simulations assume that such measurements are acquired every two days in the four monitoring columns marked in Figure 2 as “C data.” Thus, the inverse problem includes a total of 219 C data points.

Levenberg-Marquardt algorithm for parameter estimation

Our parameter estimation uses the Levenberg-Marquardt modification of the Gauss-Newton algorithm (Levenberg, 1944; Marquardt, 1963). This optimization method is based on minimizing the quadratic approximation of the regularized objective function:
Our implementation follows that of Finsterle and Kowalsky (2011), where in the first term zo=(zo,Czo,ERT) is a stacked vector of total size N, combining hydrologic (C data) and geophysical (ERT data) observations. Furthermore, Czz is the a priori covariance matrix that is a diagonal N×N matrix containing the observation errors, m is the vector of M model parameters, to be regularized by λ and W, as will be detailed shortly.
Minimizing equation 11 involves constructing the N×M Jacobian matrix
by repeated calculation of the composite vector of forward-modeling responses, z(m)=(zC(m)zERT(m)). For a given model parameter j, we obtain the corresponding matrix column (z1/mj,,zN/mj)T of J through perturbation of mj by a small quantity Δmj, thus, approximating each element through zi/mj(zi(mj+Δmj)zi(mj))/Δmj. Let us assume, permeability is our parameter of interest, so mj=kj. To quantify the effect of the perturbation kj+Δkj on the flow state, we recalculate the hydrologic forward-modeling problem given in equations 17. As it alters the Darcy filtration velocity q (equation 7), the perturbation of kj effects the distribution of solute concentration c (equations 8 and 9), which in turn will lead to perturbed electric fields E in the coupled geophysical forward-modeling problem (equation 10). The entries of the geophysical data portion of J thus take the form ziERT/mj(Ei(kj+Δkj)Ei(kj))/Δkj, where Ei is the one electric-field datum of our ERT survey array.

In equation 11, the regularization is provided by the second term, where λ is the damping parameter and W is a diagonal M×M weighting matrix that, for simplicity, we take to be the identity matrix. In principle, the regularized inverse modeling attempts to avoid high contrasts between adjacent model parameters, under the assumption that high contrasts are geologically less reasonable. In practice, we use differential smoothing by splitting up the smoothing term in equation 11 into one for horizontally and one for vertically connected model parameters. Thus, the smoothing term becomes λhmTWhm+λvmTWvm. Differential smoothing allows for the fact that horizontally stratified geologic media often exhibit larger vertical property variations, which can be accommodated by setting λh>λv.

Our numerical examples will involve a mapping of material properties between two different grids; that is, the vector m involves a distribution of parameters on a coarser (rectangular) grid than the grid underlying the actual flow-process forward modeling. The primary reason for such model scaling is typically for computational efficiency (e.g., Salazar and Piamo, 2007). For the interested reader, Appendix  B provides more details about the property-mapping procedure used here.

Excursion: Deterministic versus stochastic parameter estimation

The minimization of the objective function, as given in equation 11, defines a class of parameter estimation methods referred to as deterministic. Its counterpart is the group of stochastic inversion methods. They complement each other in that deterministic methods have the primary purpose of determining model parameters that produce a single best-fitting model. Stochastic inversions, on the other hand, produce many possible or plausible solutions, with the primary purpose of obtaining uncertainties associated with the inversion process. Most of stochastic inversion approaches operate on the basis of Bayesian statistics, in which Markov Chain Monte Carlo sampling has proven to be robust. Good overviews of this vast field including many further readings are given by Gelhar (1986), Vrugt et al. (2009), Rubin and Hubbard (2005b), Scholer et al. (2011), and Binley et al. (2015).

The increasing tendency of fluid-flow modeling toward forecasting has given importance to uncertainty quantification and thus stochastic approaches. Moreover, the petrophysical link in the joint inversion framework presented here introduces an additional layer of uncertainty, together with a higher degree of nonlinearity. Thus, stochastic methods are naturally suited because of their greater exploration of the model space and their avoidance of potential problems due to difference calculations in the presence of nonlinear functions. Most coupled hydrogeophysical works discussed above build on the stochastic approach, whereas a few use deterministic methodologies (e.g., Kowalsky et al., 2011; Vilhelmsen et al., 2014), or a hybrid (e.g. Finsterle and Kowalsky, 2008; Busch et al., 2013).

Computationally though, the stochastic method is typically orders of magnitude more expensive. Hydrogeophysical inversion applications have thus been limited to models described by only a few parameters. Because estimating complex hydrologic property distributions in a large-scale imaging style will remain a challenge to the kind of process imaging considered here, deterministic methods are likely to play an important role in the foreseeable future. Attempts at improving the efficiency of deterministic methods include the adjoint method. Under certain assumptions, it allows calculating the gradient of the cost functional in equation 11 to perform a nonlinear conjugate gradient model update (Santos et al., 2006), or the sensitivity matrix for a Gauss-Newton type of inversion (Pollock and Cirpka, 2010).

Strongly related to stochastic methods, we want to briefly allude to data assimilation, a set of statistical techniques including Bayesian statistics aimed at state forecasting by combining (assimilating) observational data and a priori knowledge (Evensen, 2009). Data assimilation methods have found their way from atmospheric forecasting into predictive surface-subsurface modeling (e.g., Houser et al., 1998), then into petroleum reservoir variable estimation (e.g., Gu and Oliver, 2007) and hydrogeophysics (e.g., Chen et al., 2012). A powerful data assimilation method, the ensemble Kalman filter (EnKF), was reviewed by Aanonsen et al. (2009) for history matching. Camporese et al. (2015) use EnKF to compare uncoupled against coupled hydrogeophysical inversion using ERT data, where they aptly use the term “imaging a process” because EnKF solves the combined parameter and state estimation problem.

Finally, recent hydrogeophysical model estimation efforts involve combinations of statistical methods referred to as machine learning (e.g., Friedel, 2016). Note that learning types of statistical algorithms have been around for a long time; only ever-increasing computational abilities are making them available for complex reservoir engineering problems. Examples are fuzzy clustering techniques (e.g., Paasche et al., 2006), artificial neural networks (e.g., Karahan and Ayvaz, 2008; Karimpouli et al., 2010), and support vector machines (Al-Anazi and Gates, 2010). See also Karpatne et al. (2018) for a general discussion on machine learning in geosciences.

Estimating permeability

We go back to our inversion demonstration, where we first estimate only the permeability field k(r). Porosity is assumed constant with ϕ(r)=0.25. Together with a constant cementation exponent m(r)=1.5, the formation factor is thus constant at F(r)=ϕm=8. We compare three end images in Figure 22. These result from inverting the, respectively, salinity data (C data), ERT data, and joined (C + ERT) data. For brevity, these inversions are called K(C), K(ERT), K(C + ERT), where K indicates the inversion for permeability. All inversions use a uniformly distributed permeability of 1011  m2 as the initial model guess. Table 1 summarizes noteworthy input and output for these three inversions in the columns under their respective names.

The preferential flow path that becomes visible in Figure 2 causes tracer accumulation due to reduced flow in the upper model region. However, Figure 2 reveals that the overall permeability structure cannot be discriminated by inverting the C data alone. The fact that the C data are depth averaged, that is, calculated from vertical column averages, explains the lack of sensitivity to property variations in the vertical direction. In contrast, the ERT electrode array geometry has the spatial resolution that is necessary to discriminate between the two major permeability regions, as can be seen from inversion K(ERT) using only ERT data (Figure 2). The joint inversion (Figure 2) yields no significant qualitative image improvement over Figure 2, whereas a slight quantitative improvement can be measured by means of a model root-mean-square (rms) error applied to the permeability results (as will be further explained below). Line 15 in Table 1 shows that this rms decreases from 3.75 for K(ERT) to 2.88 for K(C + ERT).

Viewing this example in a process-imaging context, we want to verify whether the tracer spreading based on the estimated permeability model compares to the actual case. Figure 3 captures the C field at selected times of 30, 70, 110, and 161 days, where the true C field snapshots in the upper row compare against those calculated from the three inversions (one per row). We can verify qualitatively that inversions K(ERT) and K(C + ERT) are able to reproduce the main permeability structures sufficiently well because corresponding preferential flow paths lead to a good approximation of the actual C field at the selected flow times.

For this scenario, the results demonstrate that the ERT data would provide a certain forecasting capacity. This shows in the good C field comparison at 161 days, which is 51 days after the last ERT survey time at 110 days.

To quantify the inversion performance, we first calculate the rms errors for each data component in Table 1 (rows 10–13). Comparing initial to final rms values, the fits improve by roughly an order of magnitude. However, owing to the nonuniqueness problem in inversions with insufficient data, one may end up with good data fits, whereas the model outcome differs profoundly from the true geology. Hence, synthetic inverse-modeling studies, in which the true model is known, should also supply appropriate misfit measures for the estimated parameters.

Figure 3 already indicates that the permeability models obtained through the two inversions K(ERT) and K(C + ERT) lead to a reasonable C field prediction. For a more quantitative performance measure, we use an rms measure calculated between the properties of the true model mtrue and the estimated properties mest, specifically,
where δi denotes a standard deviation. Appendix  C contains more details on the implementation of this error measure. Table 1 summarizes these rms error measures for the three inversion results of Figure 2, in which the rms calculated from the initial parameter guess (line 14: initial k model) can be compared to the end result (line 15: final k model).

Estimating permeability and formation factor

Our second study investigates inverting simultaneously for heterogeneous permeability k(r) and heterogeneous formation factor F(r). This example demonstrates a possible approach of addressing petrophysical heterogeneity, that is, allowing for a nonuniform petrophysical parameter distribution in the inversion.

Parameter sensitivity studies usually precede inversion experiments to investigate the effect on the forward-modeling response z(m) with respect to input parameter changes. A quantitative measure is given by the scaled sensitivity matrix, which consists of the elements (Finsterle, 2015)
Element Sij is the partial derivative of the calculated system response zi with respect to parameter mj, scaled by the respective standard deviations (we use the symbol δ to avoid confusion with electrical conductivity σ) δzi and δmj.
We combine the scaled sensitivity coefficients of all available observations to get the relative measure
for a given parameter j(j=1,,M). In Figure 4, the measure Sj is calculated for each cell parameter of the imaging mesh for k(r) (Figure 4) and for F(r) (Figure 4). To highlight the relative sensitivity difference between these types, Sj is not weighted, so δzi=δmj=1. The comparison reveals a sensitivity difference of up to three orders of magnitude; that is, the formation factor has a far smaller influence on the flow field evolution than does permeability.

Figure 5 shows actuals and estimations of an inversion for k(r) and F(r). This model also assumes a constant homogeneous porosity of ϕ(r)=0.25; thus, the variation in F(r) is only due to the cementation exponent m(r). Table 1 summarizes its input and output under the name KF(C + ERT), where F stands for the formation factor. The underlying true permeability model is the same as shown in Figure 2. For creating the heterogeneous model for the formation factor, we used the same unconditional sequential Gaussian simulation parameters as for permeability, yet inverting the outcome; that is, high permeability corresponds to a low F, and vice versa, as shown in Figure 5 and 5. Furthermore, the relationship between k and F was perturbed by spatial random noise Δk with a 5% maximum amplitude, specifically, F1/log10(k+Δk).

Qualitatively comparing the inversion outcome for k(r) and F(r) (Figure 5 and 5) against the true case (Figure 5 and 5) shows that the main model structures can be reproduced. However, compared to the case with known and constant F, a slightly deteriorated permeability image results leading to a deteriorated reproduction of the C field as shown in Figure 6.

Estimating permeability with the smoothing regularization term inactive

Finally, Figure 7 addresses the regularization aspect mentioned above under the comparative discussion of coupled and uncoupled inversion schemes. To recall, we want to clarify the fact that the quasiregularization imposed on the geophysical property evolution through the physics of the fluid-injection process does not extend to the permeability parameter grid in this framework of coupled processes.

To demonstrate the necessity for regularization, we repeat the joint inversion for permeability, which was called K(C + ERT). Now, the regularization parameter λ (equation 11) is set to a very small value (λ=104), thus essentially deactivating the spatial smoothing effect on the permeability parameter grid. Figure 7 compares the original result (this is the same as shown in Figure 2), against the nonregularized outcome (Figure 7), showing that the general permeability structure is still visible; however, it is more fragmented. The corresponding C field predictions in Figure 8 (Figure 88: true, Figure 88: calculated from the inversion result) also reveal the deteriorating effect.

The last column in Table 1 further quantifies the effect on the data and model rms measures. Compared to its regularized counterpart, all rms measures reveal less of an agreement with the true case. These results indicate the stabilizing nature of the smoothing term when inverting for a heterogeneous permeability distribution.

Given their more economic nature and larger volume coverage over in situ borehole observations, remote-sensing geophysical methods are likely to become increasingly important in understanding and monitoring anthropogenic subsurface alterations due to fluid injection. This didactic paper attempts to broaden the view of how to use geophysical data to image the fluid-injection process. Conventional geophysical inverse modeling does not include the forward modeling of the fluid injection and associated subsurface alterations. In contrast, the process-imaging approach put forth here inverts for the underlying process-driving (hydrologic) properties. We have demonstrated how estimates of permeability k(r) and formation factor F(r) can be obtained from a coupled inversion framework exploiting the information content of time-lapse geophysical measurements.

Estimating model parameters relevant to the flow process starts with properly forward modeling the various interactions of coupled physical systems. Owing to complex and yet to be fully understood process interactions and corresponding rock-physics relations, inverse modeling for the range of parameters involved is a highly nontrivial task. Under several simplifying assumptions, we show that time-lapse ERT data alone are sufficient for properly reproducing a salinity-injection scenario. Specifically, it was not necessary to have additional borehole measurements of salinity (what we called C data) to obtain an accurate estimate of how the salinity plume evolves through time and space.

This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division under Contract DE-AC02-05CH11231. Parts of the modeling studies were supported through a laboratory-directed research and development (LDRD) grant by Lawrence Berkeley National Lab for MC. We are grateful to John Etgen, Rune Mittet, and three anonymous reviewers for improving this work.

Data associated with this research are available and can be obtained by contacting the corresponding author.


To address the importance of diffusion to the evolution of the salinity plume, we return to equation 3 for the dispersion tensor D. We note that the first term on the right side corresponds to solute diffusion, whereas the second and third terms are due to pore-scale advection processes (those occurring within the voxels of the modeling grid) that correspond to perturbations from the average advection throughout a voxel of porous material. One can quantify the ratio of advection to diffusion by means of a local Péclet number Pe. Essentially, Pe reflects the ratio of the second to first terms in the dispersion tensor (equation 3), that is,
where γld and d is a characteristic grain diameter. Now for the ratio of the mean advection to diffusion at the macroscopic scale L of the plume (say meters or larger), we obtain a macroscopic-scale definition of the Péclet number Pe. The latter is the ratio of the advective derivative (the second term on the left side of the solute balance equation 1) to the first term on the right side, using only the diffusion contribution to the dispersion tensor, resulting in
Then, the ratio of the two Péclet numbers,
defines the importance of local solute dispersion, as controlled by length d, in comparison to macroscopic plume-scale dispersion, as controlled by length L, through the heterogeneity in the Darcy flow q in the advective derivative of the solute balance. Our numerical examples assume that the local dispersion associated with dispersivities γl and γt in the dispersion tensor are negligible relative to the macroscopic dispersion, that is, Pe/Pe1. We will keep the diffusion term Dm/F in the dispersion tensor to allow for those portions of the modeling domain that have very slow flow and, therefore, an important diffusion contribution. However, through most of the domain, we will always have Pe1. Thus, diffusion, and therefore the formation factor, will have minimal influence on the plume evolution.


Equations 18 describe a system where hydrologic state changes result from the increased flux of dissolved ions caused by electrolyte injection. Forward modeling of the evolving concentration field c(r,t) involves a (hydrologic) finite-volume mesh, ΩH, with 5098 elements, with Figure 2 indicating its mesh fidelity. The inversion domain comprises a subset of 3627 elements. To reduce such an overly fine parameter grid, the inversion parameters involve a coarser imaging mesh, ΩHi, with M=328 rectangular cells. Each cell parameter encompasses an average of 11 ΩH elements. The discretization difference between ΩH and ΩHi can be seen by comparing Figure 2 against the coarse-mesh images in Figure 22.

Each coarse-mesh cell holds a voxel value of the properties that are to be estimated, k(r) and F(r). After a model update of the iterative inversion process, these properties are mapped from ΩHi to ΩH. The changes in σf, discretized on ΩH, infer changes in the bulk rock electrical conductivity σR, calculated from equation 9. The geophysical property σR resides on another mesh ΩG, so the petrophysical transform involves a mapping from ΩH to ΩG. We use a simple nearest-neighbor scheme for this operation. All forward modeling and mesh mapping operations are modularized components of the MPiTOUGH2 software (Commer et al., 2014).

Finally, forward modeling of the ERT method is carried out on ΩG. The corresponding finite-difference Poisson-type equation and its 3D solver are described in detail by Commer et al. (2011). This ERT solver mimics a quasi-2D case by using a very coarse finite-difference grid perpendicular to the flow direction (y-axis) compared to the (x) axis of the flow direction. Preceding tests confirmed the accuracy of this approach by comparing against the 2D simulator used by Kemna et al. (2002).


When applying the rms norm of equation 13, different options exist for treating the scale differences between the true model grid ΩH (Figure 2) and the coarser grid of inversion parameters ΩHi (Figure 22). We choose a cubic spline interpolation to map properties from the coarse inversion parameter grid to the finer topology of the true model grid. The latter defines the summation in equation 13, amounting to M=3627 fine-grid differences.

Using equation 13 for our permeability outcome, we get the rms between the true (ktrue) and the estimated (kest) logarithmic (base-10) permeability field,
Table 1 summarizes these rms errors for the three inversion results of Figure 2, where the starting model rms (initial K model) can be compared to the end result (final K model). These calculations assume a constant standard deviation of δ=0.12, which is 1% of 1 darcy in log10-space.
Given that permeability has the dominating influence on preferential fluid flow pathways, the next step is to quantitatively assess the forecasting accuracy of the estimated permeability fields. The step involves forward modeling the brine injection process over the whole flow time range of 161 days. For all inversion results, we interpolate the estimated k and F fields (Figures 22, 5, 5, 7, and 7) from the coarse-grid (ΩHi) onto the fine-scale grid (ΩH) used for true-data calculation (Figure 2). Performing an outer summation over equation 13, where the variable m is now replaced by the concentration C, extends our C field imaging quality measure over the whole fluid-injection period (161 days),
Again, δ is chosen to be constant and is 1% of C=1, the maximum concentration, thus, δ=0.01. Table 1 lists all rms values calculated in equation C-2, referred to as “C field model rms,” where each initial rms (calculated from the starting model guess) can be compared against its final counterpart.
Freely available online through the SEG open-access option.