**Research article**
07 Nov 2019

**Research article** | 07 Nov 2019

# An extended trajectory-mechanics approach for calculating the path of a pressure transient: travel-time tomography

Donald W. Vasco Joseph Doetsch and Ralf Brauchler

^{1},

^{2},

^{3}

**Donald W. Vasco et al.**Donald W. Vasco Joseph Doetsch and Ralf Brauchler

^{1},

^{2},

^{3}

^{1}Energy Geosciences Division, Building 74, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA^{2}Department of Earth Sciences, ETH Zurich, Zurich, Switzerland^{3}Waste Disposal and Hydrogeology, AF-Consult Switzerland Ltd, Täfernstrasse 26, 5405 Baden, Switzerland

^{1}Energy Geosciences Division, Building 74, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA^{2}Department of Earth Sciences, ETH Zurich, Zurich, Switzerland^{3}Waste Disposal and Hydrogeology, AF-Consult Switzerland Ltd, Täfernstrasse 26, 5405 Baden, Switzerland

**Correspondence**: Donald W. Vasco (dwvasco@lbl.gov)

**Correspondence**: Donald W. Vasco (dwvasco@lbl.gov)

Received: 02 May 2019 – Discussion started: 05 Jun 2019 – Revised: 09 Sep 2019 – Accepted: 30 Sep 2019 – Published: 07 Nov 2019

The application of a technique from quantum dynamics to the governing equation for hydraulic head leads to a trajectory-based solution that is valid for a general porous medium. The semi-analytic expressions for propagation path and velocity of a change in hydraulic head form the basis of a travel-time tomographic imaging algorithm. An application of the imaging algorithm to synthetic arrival times reveals that a cross-well inversion based upon the extended trajectories correctly reproduces the magnitude of a reference model, improving upon an existing asymptotic approach. An inversion of hydraulic head arrival times from cross-well slug tests at the Widen field site in northern Switzerland captures a general decrease in permeability with depth, which is in agreement with previous studies, but also indicates the presence of a high-permeability feature in the upper portion of the cross-well plane.

Understanding the spatial variation in subsurface flow properties is important for many applications, such as groundwater extraction and storage, hydrocarbon production, geothermal energy generation, and waste water disposal. Advanced production processes like hydraulic fracturing require the development of high-resolution reservoir models necessary to capture the influence of the fractures (Zhang et al., 2014; Fujita et al., 2015). Often there are very few observations with which to infer such properties: typically measurements from a few wells intersecting a formation of interest. However, developments such as cross-well transient pressure testing (Hsieh et al., 1985; Paillet, 1993; Karasaki et al., 2000); and hydraulic tomography (Tosaka et al., 1993; Gottlieb and Dietrich, 1995; Butler et al., 1999; Yeh and Liu, 2000; Vasco and Karasaki, 2001; Bohling et al., 2002, 2007; Brauchler et al., 2003, 2010, 2011, 2013; Zhu and Yeh, 2006; Illman et al., 2007, 2008; Fienen et al., 2008; Bohling, 2009; Cardiff et al., 2009, 2013a, b; Huang et al., 2011; Sun et al., 2013; Paradis et al., 2015, 2016), have improved the ability to resolve two- and three-dimensional variations in hydraulic properties. New techniques, including fiber-optic temperature and pressure observations, and geophysical observations sensitive to pressure changes (Yeh et al., 2008; Rucci et al., 2010; Marchesini et al., 2017), will further improve spatial and temporal coverage and generate large data sets. Finally, the joint interpretation and inversion of geophysical and hydrological data leads to better constrained imaging of flow properties (Rubin et al., 1992; Hyndman et al., 1994, 2000; Vasco et al., 2001; Vasco, 2004; Kowalsky et al., 2004; Day-Lewis et al., 2006; Brauchler et al., 2012; Lochbühler et al., 2013; Soueid Ahmed et al., 2014: Ruggeri et al., 2014; Jimenez et al., 2015; Binley et al., 2015; Linde and Doetsch, 2016).

The characterization of complicated aquifer and reservoir models using sizable data sets points to the need for robust and efficient approaches for modeling pressure propagation. To this end, there are a number of approaches that aim to reduce the computational burden and data handling requirements associated with hydraulic tomography. For example, there are methods that reduce the governing equation to a simpler form for the moments of the transient head or pressure variation (Li et al., 2005; Yin and Illman, 2009; Zhu and Yeh, 2006). There are also approaches for the analysis of sinusoidal and oscillatory pumping tests that are based upon the phase shifts and amplitude differences between observed and calculated pressure variations, using these phase shifts to infer properties between two wells (Bernabe et al., 2005; Black and Kipp, 1981; Cardiff et al., 2013b; Kuo, 1972; Rasmussen et al., 2003; Renner and Messar, 2006). Another technique relies upon a measure of the arrival time of a pressure pulse or disturbance as a basis for transient travel-time imaging or tomography (Vasco et al., 2000; Kulkarni et al., 2001; Brauchler et al., 2003, 2007, 2010, 2011, 2013; He et al., 2006; Hu et al., 2011; Vasco and Datta-Gupta, 2016). Finally, there are methods that attempt to find lower-dimensional representations of the model or of the matrices describing the forward and inverse problems. These methods include principal component analysis (Lee and Kitanidis, 2014), Karhunen–Loèève expansions (Zha et al., 2018), and reduced-order models (Liu et al., 2013).

There are at least three advantages associated with the use of travel times, an alternative to the direct treatment of the entire transient head or pressure waveforms. First, the arrival of the early onset of the transient pressure pulse can be much sooner than the time at which steady-state conditions are achieved. Thus, cross-well slug tests can be conducted rapidly, facilitating improved spatial coverage. Second, the relationship between such travel times and hydraulic diffusivity is quasi-linear and convergence to a solution is not as sensitive to the initial model as it is for the direct inversion of transient pressure waveforms (Cheng et al., 2005). Third, the interpretation and reduction of transient head or pressure waveform data can be more complicated due to the sensitivity of amplitudes to various factors such as the packer coupling, the calibration of the receiver transducers, and the conditions surrounding the borehole.

Previous trajectory-based formulations of pressure arrival-time tomography relied upon an asymptotic approach that assumes smoothly varying properties (Vasco et al., 2000; Brauchler et al., 2003, 2007; He et al., 2006; Vasco, 2008; Vasco and Datta-Gupta, 2016). This assumption is certainly violated in many commonly encountered situations, such as a layered sedimentary environment and in the presence of faults or fractures. Here we apply a newly developed trajectory-based technique for travel-time tomography that dispenses with the assumption of smoothly varying properties, enlarging its range of validity to any model that may be treated using a numerical simulator (Vasco, 2018; Vasco and Nihei, 2019). The semi-analytic approach provides insight into factors controlling the propagation of a pressure transient in a complex porous medium. As shown here, the expression for the trajectories may form the basis for efficient sensitivity computations. These sensitivities are particularly useful in inverting transient pressure propagation times and in hydraulic travel-time tomography. All of the sensitivities required for the interpretation of a pressure test can be obtained in a single numerical simulation of the test. We apply the method to cross-well hydraulic tomographic imaging, considering both synthetic and field pressure arrival times.

In this section we describe our iterative algorithm for updating an aquifer model in order to improve the fit to a set of observed arrival times. We shall only discussion the elements of the derivation of Vasco (2018), as well as a perturbation technique, which are essential for understanding the inversion procedure. The approach involves a number of steps, beginning with the equation governing the transient variation in hydraulic head, and ending with a linear system of equations to be solved for the aquifer parameters. As an overview, the major steps of the methodology are shown schematically in Fig. 1. The approach is an offshoot of trajectory-based techniques developed in quantum dynamics for the study of large chemical systems (Wyatt, 2005; Liu and Makri, 2005; Goldfarb et al., 2006; Garashchuk, 2010; Garashchuk and Vazhappilly, 2010; Garashchuk et al., 2011; Gu and Garashchuk, 2016). As shown in Vasco (2018), the trajectory-mechanics treatment leads to a set of coupled ordinary differential equations that may be solved numerically, as is done in quantum mechanics. However, one can take advantage of existing numerical simulators to compute one of the unknown vector fields, reducing the system to a single set of equations for the trajectory (Vasco, 2018). The result of this analysis is a semi-analytic expression for the path of a transient pulse. This expression, along with a perturbation technique, provides a basis for an efficient method for imaging spatial variations in hydraulic diffusivity in the subsurface – a form of travel-time tomography. We illustrate the procedure with applications to both synthetic and observed arrival times in the section that follows this description.

## 2.1 Governing equation and trajectory calculations

We begin with the equation governing the evolution of a
transient variation in hydraulic head *h*(**x**,*t*) (*L*) as a function of
space **x** and time *t*,
adopting the form of the governing equation presented in de Marsily (1986, p. 109):

where **K** is the hydraulic conductivity (*L*∕*T*), a symmetric tensor,
and *ζ* is the specific storage coefficient with dimensions of the inverse length (1∕*L*).
The specific storage coefficient depends upon
the total porosity of the medium, the isothermal compressibility of the liquid,
the compressibility of the solid constituents,
and the compressibility of the porous matrix,
as discussed in de Marsily (1986, p. 109).

From this point on, we shall assume that the hydraulic head has been normalized by dividing
both sides of Eq. (1) by a constant reference head value *h*_{0} (*L*).
We will still use the variable *h*(**x**,*t*) for the normalized head, which is now unitless.
An expression for the trajectory associated with the propagation of a transient
fluid front follows from substituting the exponential representation

into the governing Eq. (1) for hydraulic head.
Because we can choose the reference location such that the hydraulic
head is always positive, Eq. (2) is well defined and can always be solved for *S*.
Upon substituting for *h*(**x**,*t*) in Eq. (1),
the resulting equation for *S*(**x**,*t*), known as the phase, may be written as

The vector ** p** (1∕

*L*) is the spatial gradient of the phase

and ** v** (

*L*∕

*T*) is a velocity vector given by

Note that Eq. (3) has the form of a traveling front with a
velocity that depends upon the vector ** p** and the
medium properties

*ζ*and

**K**. As shown in Vasco (2018), the partial differential Eq. (3) is equivalent to the system of ordinary differential equations

One can solve the two ordinary differential equations
for the trajectory **x** and the vector ** p**
(Cash and Carp, 1990; Press et al., 1992; Wyatt, 2005).
An alternative approach is to use a reservoir simulator to calculate

*h*(

**x**,

*t*) and then use Eqs. (2) and (4) to determine

**from the hydraulic head**

*p*
Substituting for ** p** in Eq. (6) gives
an expression for the trajectory in terms of

*h*(

**x**,

*t*)

We use the TOUGH2 numerical simulator (Pruess et al., 1999) to calculate the pressure and head changes and then use Eq. (9) to find the trajectories.

## 2.2 Semi-analytic sensitivities

A primary application of the trajectories described above will be to estimate flow properties between boreholes via hydraulic tomographic imaging. In this procedure a series of pumping tests are conducted in isolated segments of one borehole. During each test a rapid injection is used to generate a transient fluid pressure pulse that propagates to pressure sensors in an adjacent well. For an impulsive source, the time at which the peak pressure is observed in the adjacent borehole is defined as the arrival time. For the inverse problem we determine the flow properties from the arrival times observed in isolated sections of the monitoring well. In order to solve the inverse problem we must relate the travel time of the pressure pulse to the hydraulic properties of the medium.

Our approach to the solution of the nonlinear inverse problem will be iterative
in nature.
That is, in order to estimate flow properties we begin with an
initial model and progressively update it, solving the forward problem
of reservoir simulation at each step.
We shall need model parameter sensitivities, which are
the partial derivatives of each observation with respect to changes in
each of the model parameters (Jacquard and Jain, 1965), for every iterative update.
We will be interested in transient pressure arrival times that are
defined as the time at which the peak of a pressure pulse is
observed at a measurement point.
Expression (Eq. 9) forms the basis for our sensitivity estimates.
The only nonzero component of the velocity vector ** v** is
along the trajectory

**x**(

*t*), and it is given by the magnitude of the vector, denoted by

*v*. Integrating Eq. (9) along the path

**x**(

*t*) we have

where $x=\left|\mathbf{x}\right|$ is the distance along the path **x**.
One could relate perturbations in the arrival time of a pressure pulse to changes in the velocity but this will lead to a sensitivity
that varies as *v*^{−2}.
This will magnify the influence of any variations in velocity along the
trajectory, potentially leading to instabilities in the inversion.
Formulating the inverse problem in terms of the slowness, *s* (*T*∕*L*), given by

eliminates this problem and leads to

an integral relationship where the nonlinearity is contained entirely within
the definition of the path of integration.
That is, according to Eq. (9),
the path of integration **x** depends upon ** v** and, in turn,

*s*(

*x*).

Model parameter sensitivities, in this case relating small changes in
the slowness along the trajectory, *δ**s*(**x**),
to changes in the travel time of a transient pressure pulse,
follow from a perturbation argument.
Specifically, we consider a perturbation of the slowness with respect to
a background model *s*_{o}(**x**):

where *δ**s* is assumed to be small.
There is a corresponding small change, *δ**T*(**x**),
in the travel time from a source location to an observation point

Substituting the perturbed forms of *s*(**x**) and *T*(**x**),
given by Eqs. (13) and (14),
into the expression (Eq. 12) produces

where the integration is along a perturbed path
$\mathbf{x}={\mathbf{x}}_{o}+\mathit{\delta}\mathbf{x}$.
It has been shown that perturbations in the path lead to
terms that are second order in *δ**s*.
Thus, in computing the sensitivities, which are first order in *δ**s*,
we can neglect perturbations in the trajectory due to perturbations in *s*.
Therefore, we can integrate along the path calculated for the current or background model,
denoted by **x**_{o}.
Because the travel time in the background model, *T*_{o}, is the integral of the
background slowness function *s*_{o}(**x**) along the trajectory,
the initial terms on each side of Eq. (15) cancel and we are
left with

relating perturbations in the slowness, *δ**s*, along the trajectory to
perturbations in the arrival time, *δ**T*.

In order to update the model and the head or pressure field
using a numerical simulator, we shall need to map the
updated slowness estimates into the reservoir model parameters
*ζ* and **K**.
This cannot be done in a unique fashion and requires additional
information or assumptions.
Here, we will assume that the permeability tensor is isotropic,
so that it is of the form **K**=*K***I**,
where *K* is the scalar permeability and **I** is the identity matrix
with ones on the diagonal and zeros elsewhere.
If the inversion is part of a joint inversion of several data types
it might be possible to solve for *ζ* or *K* using
other information, such as geophysical observations.
In some formations, such as a clean sand,
it might be possible to relate the permeability to the porosity,
and to solve for the porosity uniquely in terms of the slowness.
Alternatively, as the porosity typically has a much smaller
range of variation than permeability does,
one might assume that the permeability dominates variations in *s*,
and hence solve for an effective permeability,
lumping both changes in *ζ* and permeability into changes in *K*.
It is evident from Eq. (9) that one has to correct
the estimates for variations in hydraulic head.
As shown below, we use the output of the numerical simulator, based upon the
current reservoir model, for this correction.

## 2.3 Comparison with existing asymptotic methods

Several trajectory-based methods for pressure arrival-time tomography (Vasco et al., 2000; Brauchler et al., 2003; He et al., 2006; Hu et al., 2011; Vasco and Datta-Gupta, p. 131, 2016) utilize a high-frequency asymptotic solution to the diffusion equation. A major assumption of such solutions is that the pressure variation is rapid in time (Virieux et al., 1994) or that the dominant frequencies in a Fourier transform of the trace are high. Equivalent results can be obtained if we assume that the medium properties are smoothly varying in comparison with the length scale associated with the propagating pressure transient or that parameters take on values in a particular range (Cohen and Lewis, 1967). In that case we can neglect the divergence term on the right-hand side of Eq. (3) and it reduces to an eikonal equation

where we have made use of Eqs. (4) and (5). There are efficient fast-marching methods for solving the eikonal equation (Podvin and Lecomte, 1991; Sethian, 1999; Osher and Fedkiw, 2003), that are applicable to modeling transient pressure propagation in high-resolution reservoir models (Zhang et al. 2014; Fujita et al., 2015). The eikonal equation is equivalent to a system of ordinary differential equations, the ray equations, defining the path of the transient pulse and the spatial variation of the phase (Courant and Hilbert, 1962).

From the high-frequency asymptotic solution and the ray equations,
Vasco et al. (2000) derived a semi-analytic expression,
in which the square root of the peak arrival time is given by the line integral
along the trajectory **x**_{eikonal} defined by the eikonal equation:

where

has units of $\sqrt{T}/L$,
**x**_{eikonal} signifies the trajectory resulting from
the solution of the eikonal Eq. (17), and *r* is the distance along the
trajectory.
Equation (18) is a nonlinear relationship between the travel time *T*_{peak}
and *φ* because the path **x**_{eikonal}
depends upon the spatial variation of *ζ* and *K*.
As in the previous subsection,
we can linearize the relationship by assuming a background
model and considering perturbations, or small changes, with
respect to the background model.
Because the perturbations in the path **x**_{eikonal} are
second order in perturbation of *φ*, we can write the
perturbed expression as

where $\mathit{\delta}\sqrt{{T}_{\mathrm{peak}}}$ is the perturbation in
the square root of the travel time and **x**_{0} is the
trajectory in the background medium.

In Vasco (2018) the limitations of the high-frequency asymptotic approach are discussed and illustrated. In particular, it is shown that for abrupt boundaries and sharp layers, the trajectories calculated using the eikonal equation bend too strongly into high-permeability regions of a half-space or layer. This leads to deviations in the trajectories from regions with high model parameter sensitivity, and the potential for errors when updating a simulation model. In the next section we will explore these limitations in the context of hydraulic tomography, using both synthetic and experimental data.

## 2.4 A linearized and iterative approach for tomographic imaging

A reservoir model is typically defined over a two- or three-dimensional grid that is used by a numerical reservoir simulator. For such a discrete model with properties defined on a grid of cells, and where one assumes constant values within each cell of the model, we can break up the path integrals Eqs. (16) and (20) into sums over all of the grid blocks intersected by the trajectories. For the integral Eq. (16) the discrete sum is given by

where *δ**s*_{i} is the perturbation of *s* in the *i*th grid block,
and *l*_{i} is the length of the trajectory **x**_{o} in that grid block.
Equation (21) constitutes a linear constraint on the
perturbations of *s* in the sampled grid blocks of the model
(those intersected by the path **x**_{0}).
By considering a number of sources and receiver pairs,
for example from a sequence of cross-well slug tests,
we arrive at a system of linear equations relating the
perturbations in *s* to perturbations in the observed arrival times.
We may write the system as a matrix equation

where *δ*** T** is a vector of travel-time residuals,

**M**is a matrix containing the path lengths in each grid block intersected by one or more trajectories, and

*δ*

**is a vector whose elements consist of the perturbations in each grid block (**

*s**δ*

*s*

_{i}or

*δ*

*φ*

_{i}for the

*i*th block).

In the iterative, linearized inversion scheme that we shall
adopt here, we start with an initial model, perhaps derived from well logs,
denoted by *ζ*_{0} and *K*_{0}, and calculate the background values of *v* or *φ*.
Depending upon the method, we either conduct a reservoir simulation and use
Eq. (9) to derive the trajectories, or
solve the eikonal Eq. (17) and calculate **x**_{eikonal} as in Vasco et al. (2000).
This allows us to compute the trajectories and the lengths in each grid block
and to construct the elements of the matrix **M**.
In order to update the reservoir model we need to relate the updated field *s*
to the model parameters *ζ* and *K*.
Recall that we can only resolve the ratio *ζ*∕*K* and
we cannot distinguish increases in *ζ* from
decreases in *K* or vice-versa.
In this example only the permeability varies, so we fix *ζ* and
only solve for changes in permeability.
If *ζ* also varies, we can only find an effective permeability
variation that will contain the effects of any variation in *ζ*.
The relationship follows from Eq. (9) and is given by

Because the head field *h*(**x**,*t*) is present in the integral expression, we need
to recalculate this field at each iteration.
But that calculation is already required in order to update the trajectory
**x**(*t*).

Due to errors in the data and modeling approximations, we do not
expect that the system of Eq. (22) will have an exact solution.
Thus, we seek a least squares solution in which the sum of the
squares of the residuals is minimized.
Furthermore, due to resolution and uniqueness issues,
a direct least squares solution of Eq. (22) will probably be
unstable and small errors will lead to large changes in
the estimates of *δ*** s** (Menke, 2012; Parker, 1994). Therefore, we introduce regularization or penalty terms in order to
stabilize the inverse problem.
The penalty terms seek to minimize the norm of the model update,
and to minimize the roughness of the updates, as measured by the
difference operators that mimic the second spatial derivatives
of the model, the model Laplacian (Menke, 2012).
The function that we are minimizing, Π(

*δ*

**), is the sum of the squares of the residuals, the weighted model norm and the weighted model roughness:**

*s*
where **L** is a matrix operator that mimics the second spatial
derivative of the model, *w*_{n} is the model norm weight, and
*w*_{r} is the model roughness weight.
Note that in Eq. (24) we weight all of the data uniformly.
It is possible to include a covariance matrix in order to account for
correlations between observations and variations in data quality (Tarantola, 2005).
Minimizing the quadratic function Eq. (24) with respect to the model parameters
leads to a linear system of equations for *δ*** s**:

The penalized least squares problem is solved for the perturbations, *δ*** s**,
using the least squares QR algorithm (LSQR) proposed by Paige and Saunders (1982).
With the solution in hand we then update the reservoir model.
Because the high-frequency asymptotic method only requires

*φ*, we do not need to convert back to the flow parameters

*ζ*and

*K*. Therefore, we can update the model, solve the updated eikonal equation, recompute the residuals, retrace the trajectories, calculate the sensitivities, and continue the process until the misfit to the travel times is reduced sufficiently. For the extended trajectory method we can use Eq. (23) to transform from

*s*to

*K*before updating the reservoir model and conducting another numerical simulation.

The linearized expression Eq. (25) also provides a basis for the assessment of a solution
to the inverse problem, that is, the calculation of model parameter resolution
and uncertainty (Parker, 1994; Aster et al., 2005; Menke, 2012).
Model parameter resolution estimates can be particularly useful in understanding
spatial averaging and nonuniqueness in hydrological inverse problems
(Vasco et al., 1997; Bohling, 2009; Paradis et al., 2016).
We can define model parameter resolution very simply in terms of the generalized inverse **G**^{†},
obtained from Eq. (25) by formally inverting the matrix on the left-hand side:

Hence, the parameter estimates for a given iteration, denoted by $\mathit{\delta}\widehat{\mathit{s}}$, are a linear function of the observations

Using Eq. (22) to replace *δ*** T** with

**M**

*δ*

**gives a relationship between the estimated parameters and the “true” parameters,**

*s*
where **R** is the resolution matrix (Menke, 2012), with rows that are
coefficients describing the averaging that occurs in estimating a parameter.
We can also make use of the linear relationship between the residuals and the model
parameter updates, given by Eqs. (22) and (27), to estimate an
a posteriori model parameter covariance matrix **C**_{ss} in terms of the
covariance matrix of the errors associated with the observed travel times.
In particular, if the data errors are Gaussian, characterized by the data covariance
matrix **C**_{TT}, then the model parameter covariance matrix
may be written in terms of the generalized inverse and the data covariance matrix

a consequence of the linear nature of the problem and the properties of the Gaussian distribution (Menke, 2012).

Cross-hole hydraulic travel-time tomography and cross-well slug tests are valuable approaches for imaging spatial variations in flow properties (Paillet, 1993; Yeh and Liu, 2000; Vasco and Karasaki, 2001; Bohling et al., 2002; Butler et al., 2003; Brauchler et al., 2007, 2010, 2011). Such tests can resolve features between boreholes, similar to cross-well geophysical imaging, and are directly sensitive to flow properties. In this section we set up a synthetic hydraulic tomographic test, roughly based upon a field experiment at the Widen site in Switzerland. Following that, we analyze data from the actual field experiment, using them to image the spatial variations of permeability between two shallow boreholes.

## 3.1 Synthetic hydraulic tomography test case

The overall setup of the test example is shown in Fig. 2, along with the reference model.
A set of sources in each well, denoted by filled squares and open circles,
transmit transient pressure signals to various receivers located in the adjacent borehole.
The reference distribution, a three-dimensional permeability model with a
dominantly vertical variation in properties, was generated stochastically.
That is, a uniform number generator was used to derive permeability multipliers
between 1 and 12 for each layer in the model.
A uniform random variation of 50 % was introduced within each layer
and this variation was smoothed using a three-point moving window.
The model extends an additional 5 m in the *x*, *y*, and *z* directions,
beyond the boundaries of the plane defined by the cross-well survey.

The reservoir simulator TOUGH2 (Pruess et al., 1999) was
used to model the complete set of cross-well slug tests that comprised the full synthetic experiment.
The computations were conducted using a three-dimensional mesh with
constant pressure boundary conditions, simulating a 300 s transient pressure
test for each source.
This interval provided enough time for any head variation to propagate from a source
to the receivers due to the high background permeability of $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{10}}$ m^{2}.
The large background permeability allowed us to match the rapid pulse propagation
between the boreholes that was observed during the actual Widen field experiment described below.
The initial conditions where a constant pressure of 0.616 MPa and a uniform temperature of 20 ^{∘}C.
The source-time function was defined by a jump in flow rate followed by an exponentially decreasing
rate.
The transient arrivals were defined as the time at which the rate of change in the pressure or head
reached a maximum value.
A set of synthetic arrival times were calculated using TOUGH2 simulations
and then used as a test data set for the imaging algorithm described above.
Uniform random deviates, with maximum variations of 5 % of the arrival time,
were generated using a pseudo-random number generator and added to the TOUGH2 calculated travel times.

In order to image the permeability variations between the boreholes we conducted a series of linearized inversion steps,
where we solve the system of Eq. (25) at each step.
The starting model is a uniform half-space with a permeability of $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{10}}$ m^{2}.
The model extends from 0.0 to 15.0 m laterally and from 0.0 to 15.0 m in the vertical direction.
We represent the cross-well area using a 33 by 33 grid of cells with a block size of 0.45 m and
embed this into a 15 m (in the *z* direction) thick three-dimensional model.
Each of the injection events was simulated for 300 s, even though the pressure transient
arrived at the observation points just a few seconds after the beginning of the test.
The pressure field from the simulation was used to compute the trajectories, using
the expression (Eq. 9) for the tangent vector,
and integrating it to construct the entire path.
The arrival times were calculated using both the eikonal Eq. (17)
and by post-processing the simulation results to estimate the arrival time
of the propagating transient as it reached each observation point.
The linearized iterative algorithm, where Eq. (25) is solved at each step,
was applied using both the eikonal equation and the extended trajectory approach to compute the
sensitivities in **M**.

The regularization weightings for each approach,
*w*_{n} and *w*_{r} in Eq. (25),
were estimated by trial and error.
In particular, a series of inversions were conducted for various values of *w*_{n} and
*w*_{r} and a balance was struck between satisfying the data and minimizing the
model norm and roughness.
For the eikonal-based inversion the misfit was calculated using travel times
from the eikonal equation.
For the new approach based upon the extended trajectories the travel-time misfit
was calculated using the pressures from the TOUGH2 numerical simulator.
The norm and roughness weights for the iterative eikonal inversion were
*w*_{r}=0.15 and *w*_{n}=0.15.
For the inversion utilizing the extended trajectories we set *w*_{r} and *w*_{n}
equal to 0.1 and 0.5, respectively.

At each iteration we solve for a permeability multiplier, a factor that is multiplied by the background permeability of the uniform starting model to get the estimated permeability. A total of 10 iterations for the eikonal-based algorithm took 6 s, whereas 10 iterations for the extended trajectory approach took 129 min, illustrating the computational advantage provided by an inversion approach based upon the eikonal equation. In Fig. 3 we plot the misfit reduction as a function of the number of iterations for both the high-frequency inversion algorithm (eikonal) and an inversion based upon the extended trajectories computed using Eq. (9). There is a large initial error reduction for both the inversion based upon the eikonal paths and the inversion utilizing the extended trajectories. However, as we continue updating the model and the size of the anomalies increase and the model becomes rougher, the error reduction for the two approaches diverge, and the eikonal-based updates no longer improve the fit when the reservoir simulator is used to calculate the arrival times. Note that the iterations do reduce the error calculated using the eikonal equation and the updated model, pointing to the differences between travel-time predictions made using a high-frequency asymptotic approach and using the pressure equations. This highlights the fact that the eikonal equation becomes less accurate as the model starts to violate the assumptions of a smoothly varying medium, an aspect supported by the results of Vasco (2018). The misfit reduction associated with an iterative inversion algorithm utilizing the extended trajectories is also shown in Fig. 3. In this case, the reduction is essentially monotonic and the final error is much less than that of the eikonal-based approach. The number of iterations required to attain convergence depends upon several factors. Two important elements are how close the initial model is to the final solution in model space and the level of errors in the observations that are being fit, including modeling errors. For the synthetic case considered here the level of random noise in the simulated arrival times is only 5 %. However, the modeling error becomes an issue when the asymptotic approach is no longer valid or because we assume that the permeability outside of the cross-well plane is uniform.

The final updated high-frequency solution, plotted in Fig. 4, contains higher permeabilities between about 5.5 and 7.0 m. However, the amplitude of the permeability multiplier is less than that of the reference model (Fig. 2). Furthermore, the amplitude of the high-permeability feature at around 9.0 m is underestimated, perhaps due to its narrow width of less than 1 m. The iterative inversion based upon the extended trajectories does image the two higher-permeability zones seen in the reference model (Fig. 2). The estimated amplitude of the features appears to be closer to those of the reference model but it does overestimate the permeability of the lower feature and underestimates the permeability of the upper zone.

A better idea of the differences in the magnitude of the two solutions is conveyed in Fig. 5,
where we plot the depth variation of the reference, eikonal-based, and extended trajectory-based
models.
That is, we display the depth variation of the average of the two models,
along with the upper and lower permeability multiplier values obtained in each depth interval.
It is evident that the solution provided by a conventional imaging algorithm which uses the
eikonal equation displays permeability changes with depth that are much smoother than the
reference model.
The extended approach does contain *K* multipliers that are similar in size to the those of the reference model.
Note that the exact locations of the high-permeability features do vary in depth, deviating
somewhat from the reference model shown in the leftmost panel.
This may be due to the wide (roughly 0.5 m) spacing of the source and receivers, and the
low spatial resolution of pressure data in general (Vasco et al., 1997).

Figure 6 provides more information regarding the misfit reductions for the inversions based upon the eikonal equation paths and the extended paths. It displays the calculated travel times plotted against travel times calculated using the reference model shown in Fig. 2. Both the initial travel times, calculated using the homogeneous background model used to start the inversions, and the final travel times based upon the models obtained at the conclusion of the algorithms, are shown in the plots. The initial travel-time estimates are all larger than the actual values calculated using the reference model. This is to be expected because the largest anomalies are the approximately order-of-magnitude increases associated with the upper and lower high-permeability layers in the model. The high-permeability channels promote rapid pressure propagation between the boreholes. The eikonal equation-based algorithm does reduce the average of the calculated travel times but does not lead to good fits. The inversion based upon the extended trajectories produces relatively good fits to the reference travel times.

An important aspect of the inverse problem is an assessment of the resulting model parameters and
estimates of their reliability.
As noted in Sect. 2, the calculation of two key components of the model assessment –
model parameter resolution and model parameter covariance or uncertainty – follow from the
generalized inverse **G**^{†} given by Eq. (26).
As noted in the discussion surrounding Eq. (21),
the coefficients required for the construction of the sensitivity matrix are provided by trajectory-based,
semi-analytic quantities, in particular by the ray lengths of the trajectories through each grid block
of the model.
Equations (28) and (29) provide the model parameter resolution and model parameter
covariance, respectively.
In Fig. 7 we plot the diagonal elements of the resolution matrix, and
the standard errors associated with the estimates are shown in Fig. 4.
The diagonal elements of the resolution matrix are plotted in the grid blocks to which they correspond.
Because the rows of the resolution matrix are averaging coefficients that are normalized
to have unit magnitude, if the diagonal element approaches one, the other elements approach zero.
Therefore, diagonal values near one signify well resolved parameters for those cells
and little lateral averaging between nearby grid blocks.
In general, the resolution for the test inversion is quite good and most parameters are
well determined.
Near the upper right-hand corner the resolution does approach zero due to the lack of sampling in that region.
The resolution is highest in the central region, away from the edges of the model,
due to the crossing of trajectories in those areas.
The estimated model parameter standard errors, also plotted in Fig. 7, have a very different
distribution, with larger values at the edges of the cross-well region where there are fewer
crossing rays.
We have scaled the uncertainties by the magnitude of the estimated model parameters in
order to plot them as percentages of the parameter estimates.
The data errors were of the order of 10 % of the magnitude of the travel-time
residuals that constitute the elements of the vector *δ*** T** in Eq. (27).
The lowest errors are in areas of high resolution, with the exception of the upper right-hand corner
where there is little ray coverage.
In general, the errors are less than 5 % of the magnitude of the estimated model parameters.

We end our treatment of the synthetic test with a discussion of some validation calculations,
in which additional sources were introduced to mimic independent pumping tests.
Two tests were simulated, with one source at the left edge of the model shown in Fig. 2, at
a height *Y* of 9.9 m, and the other with a source in the right borehole at a height of 10.3 m.
The travel times of transient pulses that propagate through the model,
computed using the TOUGH2 simulator, are shown in Fig. 2.
Following that, TOUGH2 was used to calculate propagation times through the models (shown in Fig. 4).
The respective arrival times through the two models are plotted against the travel times
for the reference model (Fig. 8).
The general trends of the residuals do agree, with increasing calculated travel times
following larger observed arrival times.
For the eikonal-based model there are notable deviations from the 45^{∘} line indicating perfect fit.
The calculated travel times are systematically larger than the reference times.
The estimates based upon the extended trajectory-based algorithm are much closer to the
reference times than the times from the eikonal approach are.

## 3.2 The Widen field experiment

The Widen field site, adjacent to the Thur River in northern Switzerland (Fig. 9), has been the subject of numerous geophysical and hydrological studies (Lochbühler et al., 2013). The primary goal of the work at the Widen site is to understand the hydrologic, ecologic, and biochemical effects of river restoration. The geophysical and hydrological experiments focused upon a sandy gravel aquifer that is in contact with an unrestored section of the river (Doetsch et al., 2010). The area was penetrated by a number of boreholes and is relatively well characterized. Borehole cores revealed that the roughly 7 m thick sandy gravel aquifer is overlain by a silty sand layer and that it sits atop a thick impermeable clay aquitard. Early work at the site included individual and joint inversions of cross-well seismic, radar, and electrical resistance tomography for a zoned model (Doetsch et al., 2010). The model was consistent with the three-layer structure defined by the existing boreholes. This study was followed by several others, including a cross-hole ground-penetrating radar investigation (Klotzsche et al., 2010), and three-dimensional electrical resistance tomographic (ERT) imaging of river infiltration into the site (Coscia et al., 2011, 2012). The three-dimensional ERT imaging indicates that the highest flow velocities occur in the middle of the aquifer, whereas the lowest speeds are at the base of the sequence in clay and silt-rich gravels. A joint inversion of geophysical and hydrological data (Lochbühler et al., 2013) between several well pairs was used to constrain spatial variations in reservoir storage and hydraulic conductivity. That study imaged the large-scale decrease in hydraulic conductivity with depth.

Cross-well slug interference tests, as described in Brauchler et al. (2010, 2011), were conducted at the site and are discussed in Lochbühler et al. (2013). In such tests, a near-instantaneous change in hydraulic head in a packed-off section of one well generates a fluid pressure transient in the surrounding region. Pressure transducers in isolated sections of a nearby well are used to measure the pulse that propagates between the wells. Both the travel time of the pulse and its amplitude can be used to infer hydraulic properties between the wells (Vasco et al., 2000; Brauchler et al., 2007, 2011; Vasco, 2008). Cross-well interference slug tests were conducted at two well pairs at the Widen site, as described by Lochbühler et al. (2013). The wells P2, P3, and P4 are roughly in a line that parallels the Thur River at a distance of 15 m from the river bank (Lochbühler et al.; 2013), as shown in Fig. 9. For our work we will focus on the P2–P3 well pair, where P3 is the source well and P2 is the observation well, some 3.5 m to the west. The tomographic system consists of two double-packers in each well, where the extent of the isolated regions was 0.25 m and the spacing of the intervals was 0.5 m. A suite of observed pressure variations for receivers in the observation well are shown in Fig. 10. We are interested in the propagation time of the pulse, as measured by the arrival time of the peak pressure at each observation point, which is referenced to the time at which the peak pressure is obtained in the source interval.

The overall inversion methodology was discussed and illustrated above and the details will not be repeated here.
During one step of the iterative linearized algorithm we minimize the weighted sum of the squared misfit,
the model norm, and the model roughness, as given in Eq. (24).
The requisite equations are given by the conditions that the total misfit is minimized, which is by
the equations that result from setting ∇Π equal to zero,
where the gradient is taken with respect to the components of *δ*** s**.
Thus, at each iteration we solve the set of linear equations, Eq. (25), for the perturbations in

**. The misfits Π(**

*s**δ*

**), plotted as a function of the number of updating steps in the iterative inversion algorithms, are shown in Fig. 11. The eikonal equation residuals, calculated by the reservoir simulator, tend to level off after about three iterations and decrease gradually as the inversion algorithm progresses. This may reflect the fact that as the heterogeneity increases, the eikonal paths begin to deviate from the actual trajectories, as illustrated in Vasco (2018). The match to the observations is shown in Fig. 12 for both the eikonal-based inversion and the inversion based upon the extended trajectories. The error reduction of 76 % for the extended inversion, shown in Fig. 11, is generally monotonic. The error reduction for the extended solution is significantly larger than that for the eikonal-based inversion. Both algorithms improve the fit to the observed arrival times, although considerable scatter remains in the residuals (Fig. 12).**

*s*The final models produced by the two inversion algorithms are plotted in Fig. 13. Both models display generally higher permeabilities at shallower depths, with values decreasing as the lower edge of the model is approached. The anomalies are largely horizontal, suggesting a generally layered structure, which is in agreement with previous studies (Klotzsche et al., 2010; Lochbühler et al., 2013; Jimenez et al., 2016; Somogyvari et al., 2017; Kong et al., 2018). The magnitude of the permeability variations is larger in the trajectory-mechanics-based inversion and a high-permeability layer is evident in Fig. 13. These general features are observable in the upper and lower permeability bounds plotted as a function of elevation in Fig. 14. Both models display a decrease in permeability with depth, but the variations in the eikonal-based inversion are somewhat smaller than those of the extended trajectory approach.

We can compare our results to previous work by Lochbühler et al. (2013), where a joint inversion of cross-well ground-penetrating radar travel times and hydraulic tomography (travel times and amplitudes) was discussed. In Fig. 15 the spatial variations of the logarithm of hydraulic conductivity corresponding to our inversion grid are plotted to the same color scale. These results correspond to part of Fig. 4h in Lochbühler et al. (2013). In addition, we extracted the highest and lowest permeability values as a function of depth in the inversion region and the average permeability at each elevation. All results show the same general decrease of permeability with depth in the aquifer, as the clay aquitard is approached. The variations in permeability in the extended approach are of the same order as the joint inversion result. As in the synthetic case, the magnitude of the variations in the eikonal equation inversion is smaller.

As a validation effort, we left out data from the fifth source from the bottom in Fig. 13
when conducting the inversion for the permeability multipliers.
This allowed us to use the resulting models of *K* variation shown in Fig. 13 to estimate the
travel times of pressure pulses from the source at position five to the corresponding observation points.
The resulting observed and calculated travel times from this experiment can then be used to validate the model,
as indicated in Fig. 16.
There is considerable scatter in the arrival times but the overall trend is a variation
that increases in correspondence with the observed arrival times.
The largest disagreement is between an eikonal-based arrival-time estimate and the observed value, but the
overall scatter seems comparable for the eikonal and extended methods.
The largest deviations are the large predicted travel times for the arrivals observed
at around 1.35 s.
Such long travel times might be due to the significant low-permeability values near source position five.
The source–receiver distribution is rather sparse, and there may not be sufficient redundancy to conduct
an accurate validation experiment.

Lastly, we used Eqs. (28) and (29) to calculate the diagonal elements of the matrices **R**
and **C**_{ss}, respectively, in order to assess our trajectory-based solution.
The diagonal elements of the matrices are plotted in Fig. 17, in the locations of the grid blocks that
they represent.
The diagonal coefficients of the resolution matrix are close to one if the parameter can be determined
without interference from other grid block estimates (Vasco et al., 1997).
That is the *i*th row of the resolution matrix contains averaging coefficients associated with the *i*th parameter.
The row approaches a delta-function-like distribution and the diagonal element approaches the value one
when there is little averaging with other parameters.
The diagonal elements of the resolution matrix in Fig. 17, with peak values around 0.6,
indicate moderate spatial averaging in these estimates.
In particular, there is greater averaging than in the synthetic test due to the fact that we
are only using sources situated in a single well in the field case.
The spatial averaging is greatest and the resolution poorest for the grid blocks at the edges of the model,
particular at the top of the cross-well region.
Similarly, the model errors, also shown in Fig. 17, are larger than in the synthetic test – around 20 %
of the size of the model estimates.
The resolution and covariance estimates indicate that the high-permeability layer,
located in the upper portion of the model,
is moderately well constrained by the observations.
Due to sampling issues the error estimates are not reliable at the edges of the model and tend to zero
where there are few or no trajectories.
As indicated in the synthetic test, putting sources in both wells would increase the resolution and
reduce the uncertainty associated with our estimates, suggesting how we might improve our imaging in the
future.

The trajectory-mechanics approach described in Vasco (2018) and applied here is very general and can be used to model other hydrological processes such as tracer transport (Vasco et al., 2018a) and multiphase fluid flow. One advantage associated with transient pressure is the rapid propagation of a disturbance in comparison with the velocities associated with fluid transport. Thus, transient cross-well pressure testing can be conducted relatively rapidly in formations with moderate hydraulic conductivity. This is particularly true when transient pressure travel times, such as the arrival time of the peak of a pressure pulse or the peak of the time derivative of the pressure (Vasco et al., 2000) are used. For the Widen field experiment the peaks are observed in the first few seconds of the measured traces in the adjacent borehole. Another advantage of hydraulic travel-time tomography is that the relationship between the arrival times and the hydraulic conductivity or diffusivity is quasi-linear (Cheng et al., 2005). Thus, the final model resulting from an inversion of travel times is less sensitive to the initial or starting aquifer model and less likely to become trapped in a local minimum. Finally, travel-time tomography provides an element of data reduction, from an entire transient pressure waveform, to a single arrival time. This can be advantageous when dealing with many intervals from multiple boreholes, time-lapse pressure changes, or large data sets derived from geophysical observations.

We have presented two examples of hydraulic tomographic imaging, one using synthetic transient pressure arrival times and the other using data from an experiment at the Widen field site on the Thur River in northern Switzerland. We do find that an algorithm based upon the eikonal equation is significantly faster than one utilizing the extended trajectories calculated using a reservoir simulator, taking only about 10 s compared with 129 min. From the synthetic application we find that an imaging technique based upon the eikonal equation, the current method used for trajectory-based modeling, has difficulty accurately imaging large and abrupt changes in permeability. Such rapid spatial changes in flow properties are a common occurrence in geologic media, with the presence of layering and fractures, with correspondingly large variations in hydraulic conductivity. For example, in well logs it is quite common to observe thin layers with permeabilities that are orders of magnitude larger than values in the surrounding formations. Indeed, in our field case at the Widen field site we image an order of magnitude change in permeability in agreement with previous results at the site. While the eikonal equation is much faster and can recover large-scale spatial variations, it is likely to produce smoothed images of sharp features and to underestimate rapid changes in properties. Thus, the approach is useful as a rapid reconnaissance tool, as in real-time imaging, and for regions where the properties are thought to be smoothly varying. This usage is supported by that fact that both the eikonal-based and the extended trajectory-based methods share the quasi-linearity of travel-time inversion approaches (Cheng et al., 2005), and are less sensitive, in comparison with inversions based upon head magnitudes, to the initial or starting model.

For a full analysis and interpretation of field data, however, we recommend the trajectory-mechanics approach; this is due to the fact that it does not invoke assumptions about model smoothness
and is therefore more robust and accurate, yet it retains the semi-analytic
sensitivities that are characteristic of trajectory-based approaches.
The semi-analytic sensitivities are computed after a single simulation, using either
numerical methods to solve the coupled system for **x** and ** p** or
using a numerical simulator to determine

**directly. Even if one resorts to a numerical simulation, the semi-analytical nature of the sensitivities provide some advantages over conventional methods. The most efficient conventional method for computing numerical sensitivities is based upon adjoint methods and requires the formulation and solution of the adjoint equation along with an additional simulation to calculate the residuals. Thus, two simulations are required in order to estimate the sensitivities for a given test.**

*p*The approach that we have described is useful for imaging permeability
variations between boreholes but it does have some limitations.
The use of slug tests limits the allowable distance between wells that
may be used for imaging variations in *K*.
However, as noted in Vasco et al. (2000), one can use
a constant rate test and consider the arrival time of the steepest slope,
extending the range of the test to larger offsets between wells.
We have chosen to fix the reservoir storage and determine variations in
an effective *K*.
This assumption needs to be explored in future studies and tested under
realistic conditions.
The computation aspects of this approach are significant, requiring
full reservoir simulations for the inversion.
As noted in Vasco (2018), there are more efficient methods
that involve solving the equations for the trajectory directly, without
a reservoir simulation.
This should reduce the computation burden of the approach at the
cost of a more complicated implementation.

The pressure data from the Widen field test are available on the Zenodo archive (https://doi.org/10.5281/zenodo.1445756, Vasco et al., 2018b).

DWV devised, implemented, and executed the inversion algorithm. JD provided and helped with the data analysis, worked on the comparison with previous studies, and participated in writing the paper. RB conceived and carried out the field experiments at the Widen field site, provided the observations and data reduction, and helped to write the paper.

The authors declare that they have no conflict of interest.

Work performed at Lawrence Berkeley National Laboratory was supported by the US Department of Energy under contract number DE-AC02-05-CH11231, Office of Basic Energy Sciences of the US Department of Energy.

This research has been supported by the US Department of Energy, Office of Basic Energy Sciences (contract no. DE-AC02-05-CH11231).

This paper was edited by Gerrit H. de Rooij and reviewed by two anonymous referees.

Aster, R. C., Borchers, B., and Thurber, C. H.: Parameter Estimation and Inverse Probl., Elsevier, Amsterdam, 2005.

Bernabe, Y., Mok, U., and Evans, B.: A note on the oscillating flow method for measuring rock permeability, Int. J. Rock Mech. Min., 43, 311–316, 2005.

Binley, A., Hubbard, S. S., Huisman, J. A., Revil, A., Robinson, D. A., Singha, K., and Slater, L. D.: The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales, Water Resour. Res., 51, 3837–3866, https://doi.org/10.1002/2015WR017016, 2015.

Black, J. H. and Kipp, K. L.: Determination of hydrogeological parameters using sinusoidal pressure tests: A theoretical appraisal, Water Resour. Res., 17, 686–692, 1981.

Bohling, G. C.: Sensitivity and resolution of tomographic pumping tests in an alluvial aquifer, Water Resour. Res., 45, W02420, https://doi.org/10.1029/2008WR007249, 2009.

Bohling, G. C., Zhan, X., Butler, J. J., and Zheng, L.: Steady shape analysis of tomographic pumping tests for characterization of aquifer heterogeneities, Water Resour. Res., 38, 60–61, 2002.

Bohling, G. C., Butler, J. J., Zhan, X., and Knoll, M. D.: A field assessment of the value of steady-shape hydraulic tomography for characterization of aquifer heterogeneities, Water Resour. Res., 43, W05430, https://doi.org/10.1029/2006WR004932, 2007.

Brauchler, R., Liedl, R., and Dietrich, P.: A travel time based hydraulic tomographic approach, Water Resour. Res., 39, 1–12, 2003.

Brauchler, R., Cheng, J.-T., Dietrich, P., Everett, M., Johnson, B., Liedl, R., and Sauter, M.: An inversion strategy for hydraulic tomography: Coupling travel time and amplitude inversion, J. Hydrol., 345, 184–198, 2007.

Brauchler, R., Hu, R., Vogt, T., Al-Halbouni, D., Heinrichs, T., Ptak, T., and Sauter, M.: Cross-well slug interference tests: An effective characterization method for resolving aquifer heterogeneity, J. Hydrol., 384, 33–45, 2010.

Brauchler, R., Hu, R., Dietrich, P., and Sauter, M.: A field assessment of high-resolution aquifer characterization based on hydraulic travel time and hydraulic attenuation tomography, Water Resour. Res., 47, W03503, https://doi.org/10.1029/2010WR009635, 2011.

Brauchler, R., Doetsch, J., Dietrich, P., and Sauter, M.: Derivation of site-specific relationships between hydraulic parameters and p-wave velocities based on hydraulic and seismic tomography, Water Resour. Res., 48, W03531, https://doi.org/10.1029/2011WR010868, 2012.

Brauchler, R., Hu, R., Hu, L., Jimenez, S., Bayer, P., Dietrich, P., and Ptak, T.: Rapid field application of hydraulic tomography for resolving aquifer heterogeneity in unconsolidated sediments, Water Resour. Res., 49, 2013–2024, 2013.

Butler, J. J., McElwee, C. D., and Bohling, G. C.: Pumping tests in networks of multilevel sampling wells: motivation and methodology, Water Resour. Res., 35, 3553–3560, 1999.

Butler, J. J., Garnett, E. J., and Healey, J. M.: Analysis of slug tests in formations of high hydraulic conductivity, Groundwater, 41, 620–630, 2003.

Cardiff, M., Barrash, W., Kitanidis, P. D., Malama, B., Revil, A., Straface, S., and Rizzo, E.: A potential-based inversion of unconfined steady-state hydraulic tomography, Ground Water, 47, 259–270, https://doi.org/10.1111/j.1745-6584.2008.00541.x, 2009.

Cardiff, M., Barrash, W., and Kitanidis, P. K.: Hydraulic conductivity imaging from 3-D transient hydraulic tomography at several pumping/observation densities, Water Resour. Res., 49, 7311–7326, https://doi.org/10.1002/wrcr.20519, 2013a.

Cardiff, M., Bakhos, T., Kitanidis, P. K., and Barrash, W.: Aquifer heterogeneity characterization with oscillatory pumping: Sensitivity analysis and imaging potential, Water Resour. Res., 49, 5395–5410, https://doi.org/10.1002/wrcr.20356, 2013b.

Cash, J. R. and Carp, A. H.: A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, ACM Transactions on Mathematical Software, 16, 201–222, 1990.

Cheng, H., He, Z., and Datta-Gupta, A.: A comparison of travel-time and amplitude matching for field-scale production data integration: sensitivity, non-linearity, and practical implications, SPE Journal, 10, 75–90, 2005.

Cohen, J. K. and Lewis, R. M.: A ray method for the asymptotic solution of the diffusion equation, J. I. Math. Appl., 3, 266–290, 1967.

Coscia, I., Greenhalgh, S. A., Linde, N., Doetsch, J., Marescot, L., Gunther, T., Vogt, T., and Green, A. G.: 3D crosshole ERT for aquifer characterization and monitoring of infiltrating river water, Geophysics, 76, G49–G59, https://doi.org/10.1190/1.3553003, 2011.

Coscia, I., Linde, N., Greenhalgh, S. A., Vogt, T., and Green, A.: Estimating traveltimes and groundwater flow patterns using 3D time-lapse crosshole ERT imaging of electrical resistivity fluctuations induced by infiltrating river water, Geophysics, 77, E239–E250, https://doi.org/10.1190/GEO2011-0328.1, 2012.

Courant, R. and Hilbert, D.: Methods of Mathematical Physics, John Wiley and Sons, New York, 1962.

Day-Lewis, F. D., Lane Jr., J. W., and Gorelick, S. M.: Combined interpretation of radar, hydraulic, and tracer data from a fractured-rock aquifer near Mirror Lake, New Hampshire, USA, Hydrogeol. J., 14, 1–14, https://doi.org/10.1007/s10040-004-0372-y, 2006.

de Marsily, G.: Quantitative Hydrogeology, Academic Press, San Diego, 1986.

Doetsch, J., Linde, N., Coscia, I., Greenhalgh, S. A., and Green, A. G.: Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data, Geophysics, 75, G53–G64, https://doi.org/10.1190/1.3496476, 2010.

Fienen, M. N., Clemo, T., and Kitanidis, P. K.: An interactive Bayesian geostatistical inverse protocol for hydraulic tomography, Water Resour. Res., 44, 1–19, https://doi.org/10.1029/2007WR006730, 2008.

Fujita, Y., Datta-Gupta, A., and King, M. J.: A comprehensive reservoir simulator for unconventional reservoirs based on the Fast Marching Method and diffusive time of flight, SPE Journal, 21, https://doi.org/10.2118/173269-PA, 2015.

Garashchuk, S.: Quantum trajectory dynamics in imaginary time with the momentum-dependent quantum potential, J. Chem. Phys., 132, 014112, https://doi.org/10.1063/1.3289728, 2010.

Garashchuk, S. and Vazhappilly, T.: Multidimensional quantum trajectory dynamics in imaginary time with approximate quantum potential, J. Phys. Chem., 114, 20595–20602, https://doi.org/10.1021/jp1050244, 2010.

Garashchuk, S., Mazzuca, J., and Vazhappilly, T.: Efficient quantum trajectory representation of wavefunctions evolving in imaginary time, J. Chem. Phys., 135, 034104, https://doi.org/10.1063/1.3610165, 2011.

Goldfarb, Y., Degani, I., and Tannor, D. J.: Bohmian mechanics with complex action: A new trajectory-based formulation for quantum mechanics, J. Chem. Phys., 125, 1–4, 2006.

Gottlieb, J. and Dietrich, P.: Identification of the permeability distribution in soil by hydraulic tomography, Inverse Probl., 11, 353–360, https://doi.org/10.1088/0266-5611/11/2/005, 1995.

Gu, B. and Garashchuk, S.: Quantum dynamics with Gaussian bases defined by quantum trajectories, J. Phys. Chem., 120, 3023–3031, https://doi.org/10.1021/acs.jpca.5b10029, 2016.

He, Z., Datta-Gupta, A., and Vasco, D. W.: Rapid inverse modeling of pressure interference tests using trajectory-based traveltime and amplitude sensitivities, Water Resour. Res., 42, 1–15, 2006.

Hsieh, P. A., Neuman, S. P., Stiles, G. K., and Simpson, E. S.: Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media, 2 Methodology and application to fractured rocks, Water Resour. Res., 21, 1667–1676, 1985.

Huang, S.-Y., Wen, J.-C., Yeh, T.-C. J., Lu, W., Juan, H.-L., Tseng, C.-M., Lee, J.-H., and Chang, K.-C.: Robustness of joint interpretation of sequential pumping tests: numerical and field experiments, Water Resour. Res., 47, W10530, https://doi.org/10.1029/2011WR010698, 2011.

Hu, R., Brauchler, R., Herold, M., and Bayer, P.: Hydraulic tomography analog outcrop study: Combining travel time and steady shape inversion, J. Hydrol., 409, 350–362, 2011.

Hyndman, D. W., Harris, J. M., and Gorelick, S. M.: Coupled seismic and tracer test inversion for aquifer property characterization, Water Resour. Res., 30, 1965–1977, 1994.

Hyndman, D. W., Harris, J. M., and Gorelick, S. M.: Inferrinng the relation between seismic slowness and hydraulic conductivity in heterogeneous aquifers, Water Resour. Res., 36, 2121–2132, 2000.

Illman, W. A., Liu, X., and Craig, A.: Steady-state hydraulic tomography in a laboratory aquifer with deterministic heterogeneity: multimethod and and multiscale validation of hydraulic conductivity tomograms, J. Hydrol., 341, 222–234, https://doi.org/10.1016/j.jhydrol.2007.05.011, 2007.

Illman, W. A., Craig, A. J., and Liu, X.: Practical issues in imaging hydraulic conductivity through hydraulic tomography, Ground Water, 46, 120–132, https://doi.org/10.1111/j.1745-6584.2007.00374.x, 2008.

Jacquard, P. and Jain, C.: Permeability distribution from field pressure data, Society of Petroleum Engineering Journal, 5, 281–294, 1965.

Jimenez, S., R. Brauchler, Hu, R., Hu, L., Schmidt, S., Ptak, T., and Bayer, P.: Prediction of solute transport in a heterogeneous aquifer utilizing hydraulic conductivity and specific storage tomograms, Water Resour. Res., 51, 5504–5520, 2015.

Jimenez, S., Mariethoz, G., Brauchler, R., and Bayer, P.: Smart pilot points using reversible-jump Markov-chain Monte Carlo, Water Resour. Res., 52, 3966–3983, https://doi.org/10.1002/2015WR017922, 2016.

Karasaki, K., Freifeld, B., Cohen, A., Grossenbacher, K., Cook, P., and Vasco, D.: A multidisciplinary fractured rock characterization study at the Raymond field site, Raymond California, J. Hydrol., 236, 17–34, 2000.

Klotzsche, A., van der Kruk, J., Meles, G. A., Doetsch, J., Maurer, H., and Linde, N.: Full-waveform inversion of cross-hole ground-penetrating radar data to characterize a gravel aquifer close to the Thur River, Switzerland, Near Surf. Geophys., 8, 635–649, https://doi.org/10.3997/1873-0604.2010054, 2010.

Kong, X.-Z., Deuber, C. A., Kittila, A., Somogyvari, M., Mikutis, G., Bayer, P., Stark, W. J., and Saar, M. O.: Tomographic reservoir imaging with DNA-labeled silica nanotracers: The first field validation, Envir. Sci. Tech., 52, 13681–13689, https://doi.org/10.1021/acs.est.8b04367, 2018.

Kowalsky, M. B., Finsterle, S., and Rubin, Y.: Estimating flow parameter distributions using ground-penetrating radar and hydrological measurements during transient flow in the vadose zone, Adv. Water Resour., 27, 583–599, https://doi.org/10.1016/j.advwatres.2004.03.003, 2004.

Kuo, C.: Determination of reservoir properties from sinusoidal and multirate flow tests in one or more wells, Society of Petroleum Engineering Journal, 12, 499–507, 1972.

Kulkarni, K. N., Datta-Gupta, A., and Vasco, D. W.: A streamline approach for integrating transient pressure data into high-resolution reservoir models, SPE Journal, 6, 273–282, 2001.

Lee, J. and Kitanidis, P. K.: Large-scale hydraulic tomography and joint inversion of head and tracer data using the Principal Component Geostatistical Approach (PCGA), Water Resour. Res., 50, 5410–5427, https://doi.org/10.1002/2014WR015483, 2014.

Li, W., Nowak, W., and Cirpka, O. A.: Geostatistical inverse modeling of transient pumping tests using temporal moments of drawdown, Water Resour. Res., 41, W08403, https://doi.org/10.1029/2004WR003874, 2005.

Linde, N. and Doetsch, J.: Joint inversion in hydrogeophysics and near-surface geophysics, in: Integrated Imaging of the Earth, edited by: Moorkamp, M., Lelievre, P. G., Linde, N., and Khan, A., 218, 119–135, AGU Geophysical Monograph Series, John Wiley and Sons Inc, Hoboken, NJ, 2016.

Liu, J. and Makri, N.: Bohm's formulation in imaginary time: estimation of energy eigenvalues, Mol. Phys., 103, 1083–1090, https://doi.org/10.1080/00268970512331339387, 2005.

Liu, X., Zhou, Q., Birkholzer, J., and Illmann, W. A.: Geostatistical reduced-order models in underdetermined inverse problems, Water Resour. Res., 49, 6587–6600, 2013.

Lochbühler, T., Doetsch, J., Brauchler, R., and Linde, N.: Structure-coupled joint inversion of geophysical and hydrological data, Geophysics, 78, ID1–ID14, https://doi.org/10.1190/GEO2012-0460.1, 2013.

Marchesini, P., Ajo-Franklin, J. B., and Daley, T. M.: In situ measurement of velocity-stress sensitivity using crosswell continuous active-source seismic monitoring, Geophysics, 82, D319–D326, https://doi.org/10.1190/GEO2017-0106.1, 2017.

Menke, W.: Geophysical Data Analysis: Discrete Inverse Theory, Academic Press, San Diego, 2012.

Osher, S. and Fedkiw, R.: Level Set Methods and Dynamic Implicit Surfaces, Springer, New York, 2003.

Paige, C. C. and Saunders, M. A.: LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Software, 8, 43–71, 1982.

Paillet, F. L.: Using borehole geophysics and cross-borehole flow testing to define connections between fracture zones in bedrock aquifers, J. Appl. Geophys., 30, 261–279, 1993.

Paradis, D., Gloaguen, E., Lefebvre, R., and Giroux, B.: Resolution analysis of tomographic slug tests head data: two-dimensional radial case, Water Resour. Res., 51, 2356–2376, https://doi.org/10.1002/2013WR014785, 2015.

Paradis, D., Lefebvre, R., Gloaguen, E., and Giroux, B.: Comparison of slug and pumping tests for hydraulic tomography experiments: a practical perspective, Environ. Earth Sci., 75, 1–13, https://doi.org/10.1007/s12665-016-5935-4, 2016.

Parker, R. L.: Geophysical Inverse Theory, Princeton University Press, Princeton, 1994.

Podvin, P. and Lecomte, I.: Finite-difference computation of traveltimes in very contrasted velocity models: A massively parallel approach and its associated tools, Geophys. J. Int., 105, 271–284, 1991.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical Recipes, Cambridge University Press, Cambridge, 1992.

Pruess, K., Oldenburg, C., and Moridis, G.: TOUGH2 User's Guide, Version 2.0, LBNL Report, 43134, Berkeley, 1999.

Rasmussen, T. C., Haborak, K. G., and Young, M. H.: Estimating aquifer hydraulic properties using sinusoidal pumping at the Savannah River Site, South Carolina, USA, Hydrogeol. J., 11, 466–482, 2003.

Renner, J. and Messar, M.: Periodic pumping tests, Geophys. J. Int., 167, 479–493, 2006.

Rubin, Y., Mavko, G., and Harris, J. M.: Mapping permeability in heterogeneous aquifers using hydrological and seismic data, Water Resour. Res., 28, 1192–1200, 1992.

Rucci, A., Vasco, D. W., and Novali, F.: Fluid pressure arrival-time tomography: Estimation and assessment in the presence of inequality constraints with an application to production at the Krechba field, Algeria, Geophysics, 75, O39–O55, https://doi.org/10.1190/1.3493504, 2010.

Ruggeri, P., Gloaguen, E., Lefebvre, R., Irving, J., and Holliger, K.: Integration of hydrological and geophysical data beyond the local scale: Application of Bayesian sequential simulation to field data from the Saint-Lambert-de-Lauzon site, Quebec, Canada, J. Hydrol., 514, 271–280, 2014.

Sethian, J. A.: Level Set Methods and Fast Marching Methods, Cambridge University Press, Cambridge, 1999.

Somogyvari, M. and Bayer, P.: Field validation of thermal tracer tomography for reconstruction of aquifer heterogeneity, Water Resour. Res., 53, 5070–5084, https://doi.org/10.1002/2017WR020543, 2017.

Soueid Ahmed, A., Jardani, A., Revil, A., and Dupont, J. P.: Hydraulic conductivity field characterization from the joint inversion of hydraulic heads and self-potential data, Water Resour. Res., 50, 3502–3522, 2014.

Sun, R., Yeh, T.-C. J., Mao, D., Jin, M., Lu, W., and Hao, Y.: A temporal sampling strategy for hydraulic tomography analysis, Water Resour. Res., 49, 3881–3896, https://doi.org/10.1002/wrcr.20337, 2013.

Tarantola, A.: Inverse Problem Theory, Society of Industrial and Applied Mathematics, Philadelphia, 2005.

Tosaka, H., Masumoto, K., and Kojima, K.: Hydropulse tomography for identifying 3-D permeability distribution in high level radioactive waste management, Proceedings of the 4th Annual International Conference of the American Society of Civil Engineers, Reston, Virgina, 995–959, 1993.

Vasco, D. W.: Estimation of flow properties using surface deformation and head data: A trajectory-based approach, Water Resour. Res., 40, W10104, https://doi.org/10.1029/2004WR003272, 2004.

Vasco, D. W.: Zeroth-order inversion of transient pressure observations, Inverse Probl., 24, 1–21, https://doi.org/10.1088/0266-5611/24/2/025013, 2008.

Vasco, D. W.: An extended trajectory mechanics approach for calculating the path of a pressure transient: Derivation and illustration, Water Resour. Res., 54, 1–19, https://doi.org/10.1002/2017WR021360, 2018.

Vasco, D. W. and Datta-Gupta, A.: Subsurface Fluid Flow and Imaging, Cambridge University Press, Cambridge, 2016.

Vasco, D. W. and Karasaki, K.: Inversion of pressure observations: an integral formulation, J. Hydrol., 253, 27–40, 2001.

Vasco, D. W. and Nihei, K.: Broad-band trajectory mechanics, Geophys. J. Int., 216, 745–759, https://doi.org/10.1093/gji/ggy435, 2019.

Vasco, D. W., Datta-Gupta, A., and Long, J. C. S.: Resolution and uncertainty in hydrological characterization, Water Resour. Res., 33, 379–397, https://doi.org/10.1029/96WR03301, 1997.

Vasco, D. W., Keers, H., and Karasaki, K.: Estimation of reservoir properties using transient pressure data: An asymptotic approach, Water Resour. Res., 36, 3447–3465, 2000.

Vasco, D. W., Karasaki, K., and Kishida, K.: A coupled inversion of pressure and surface displacement, Water Resour. Res., 37, 3071–3089, 2001.

Vasco, D. W., Pride, S. R., Zahasky, C., and Benson, S. M.: Calculating trajectories associated with solute transport in a heterogeneous medium, Water Resour. Res., 54, 1–19, https://doi.org/10.1029/2018WR023019, 2018a.

Vasco, D. W., Doetsch, J., and Brauchler, R.: Widen Field Test Pressure Data – P02 Experiment, Data set, Zenodo, https://doi.org/10.5281/zenodo.1445756, 2018b.

Virieux, J., Flores-Luna, C., and Gibert, D.: Asymptotic theory for diffusive electromagnetic imaging, Geophys. J. Int., 119, 857–868, 1994.

Wyatt, R. E.: Quantum Dynamics with Trajectories, Springer, New York, 2005.

Yeh, T.-C. J. and Liu, S.: Hydraulic tomography: development of a new aquifer test method, Water Resour. Res., 36, 2095–2105, 2000.

Yeh, T.-C. J., Lee, C. H., Hsu, K. C., and Wen, J. C.: Fusion of hydrologic and geophysical tomographic surveys, Geosci. J., 12, 159–167, 2008.

Yin, D. and Illman, W. A.: Hydraulic tomography using temporal moments of drawdown recovery data: A laboratory sandbox study, Water Resour. Res., 45, W01502, https://doi.org/10.1029/2007WR006623, 2009.

Zha, Y., Yeh, T.-C. J., Illman, W. A., Zheng, W., Zhang, Y., Sun, F., and Shi, L.: A reduced-order successive linear estimator for geostatistical inversion and its application in hydraulic tomography, Water Resour. Res., 54, 1616–1632, https://doi.org/10.1002/2017WR021884, 2018.

Zhang, Y., Bansal, N., Fujita, Y., Datta-Gupta, A., King, M. J., and Sankaran, S.: From streamlines to Fast Marching: Rapid simulation and performance assessment of shale gas reservoirs using diffusive time of flight as a spatial coordinate, SPE Journal, 21, 1–16, https://doi.org/10.2118/168997-PA, 2014.

Zhu, J. and Yeh, T.-C. J.: Analysis of hydraulic tomography using temporal moments of drawdown recovery data, Water Resour. Res., 42, W02403, https://doi.org/10.1029/2005WR004309, 2006.