Aquifer Characterization

Aquifer Characterization¹

This example illustrates how you can use COMSOL Multiphysics’ Optimization interface in combination with a PDE or physics interface to solve inverse-modeling problems (sometimes referred to as parameter estimation or history-matching problems). The modeling techniques presented here in the context of flow in an aquifer with spatially variable hydraulic conductivity are generally applicable for solving underdetermined optimization problems with COMSOL Multiphysics.

This model requires the Optimization Module.

Inverse-Modeling Background

Consider data consisting of the following components:

•

A model — for example, a PDE — of a natural or engineered system with a set of unknown coefficients or parameters

•

A set of measurements or observations with associated measurement errors

Inverse modeling is the practice of using data of this kind as input to estimate the unknown parameters of the example, which in this context is called the forward model. In particular, if the number of unknown parameters, n, is larger than the number of measurement values, m, the inverse problem is called underdetermined. This example belongs to this class of inverse problems.

When, as in this case, the forward model is a PDE and the unknown coefficient is a spatially varying random field, the number of unknown parameters is infinite. Even an ordinary finite element discretization typically gives too many unknown parameters by some orders of magnitude. Obtaining a tractable problem that can be solved numerically instead requires a separate regularization.

For underdetermined inverse problems, data fitting alone is not sufficient to single out an optimal solution; achieving this aim requires some additional criterion and an associated penalty function that ranks the solutions according how well they satisfy the criterion in question.

The objective function for an underdetermined inverse problem can thus be written as the sum of a fitness term and a penalty term:

The fitness term

(1)

measures how well the model fits with the observations. Here y denotes an m-dimensional row vector of measurement values; s is an n-dimensional row vector of parameter values; h: Rn → Rm is the forward model, which maps from parameter values to expected measurements; and R is the m-by-m covariance matrix of measurement errors. Assuming, that the measurement errors are independent and identically distributed with variance σR2, it follows that R = σR2I, where I denotes the m-by-m identity matrix.

The penalty term is relevant for problems where the number of parameters, n, exceeds the number of measurement values, m. It serves to discriminate between solutions with comparable fitness values. Through geostatistical reasoning, Kitanidis (Ref. 2) arrives at the penalty term

(2)

where Q−1 is the inverse of the spatial covariance matrix Q ≡ E[(s − Xβ)(s − Xβ)T], with E[ ] denoting the expectation value, X being an n-dimensional row vector whose elements all equal 1, and β referring to the scalar constant mean of the parameter field.

The covariance Q is a symmetric n-b-n matrix, which must be estimated from the measurements or assumed known based on prior information about the process generating s. It is typically assumed that s is a realization of a stationary and isotropic random field, such that elements of the covariance matrix are only a function of the distance between the corresponding points in space: Qij = q(|xi−xj|). The dependence of q(|xi−xj|) on distance is usually expressed as a variogram γ(xi−xj|) = q(0)−q(|xi−xj|), which is in turn assumed to have a simple functional form with only a few parameters that can be fitted to the data.

Model Definition

Forward Model

Consider a square zone in a confined 2D aquifer of side b = 100 m and with spatially variable hydraulic conductivity, Ks (m/s), which is of interest for transport modeling or other purposes. In this simplified model, assume that the region of interest is surrounded by a practically infinite domain known to have roughly constant hydraulic conductivity, Ks0. Using infinite elements, you can accurately simulate the effect of the infinite exterior domain on the region of interest.

Neglecting any differences in elevation, the hydraulic potential governing the flow through the aquifer and its surroundings can be represented by the pressure head H ≡ p/(ρf g) (m), where p is the fluid pressure (Pa), ρf (kg/m3) is the (constant) fluid density, and g (m/s2) is the acceleration due to gravity. The pressure head, or the hydraulic head, obeys Darcy’s Law

where Qs (1/s) represents sources and sinks. At the modeling domain’s boundary, on the outside of the infinite element zone, any effects of spatial inhomogeneities are negligible and you specify the condition H = 0 to complete the forward model. Note that this choice implies that the variable H refers to the change in hydraulic head from a reference state where no pumping occurs rather than the hydraulic head itself.

Inverse Model

The inverse problem is to estimate the hydraulic-conductivity field in the aquifer using experimental data in the form of hydraulic-head measurements from four dipole-pump tests. More specifically, consider eight wells (represented as points) at the corners and midpoints of the zone being characterized. While injecting fluid at one point, fluid is extracted at the same rate at the opposite point on the rim across the aquifer, and the hydraulic-head values at the six remaining points are recorded. Repeating this procedure for the other three cross-aquifer point pairs, gives a total of 4 sets of 6 observations each. The hydraulic-head measurements are assumed to be independent and have an accuracy on the order of ΔH = 1 cm , resulting in the fitness term (see Equation 1)

(3)

To make use of the experimental data, you must reduce the infinite number of degrees of freedom in the hydraulic-conductivity field to a finite number of unknown parameters. To this end, decompose the quadratic aquifer domain into a 10-by-10 grid of squares of side c = 10 m, and assume that the hydraulic conductivity takes a constant value in each square. The number of parameters characterizing the aquifer is then 100, which gives a well-defined underdetermined inverse model that can be solved quickly as an example problem.

A common assumption in the geological sciences is that spatially distributed parameters follow a geostatistical distribution defined by some parameterized variogram. This COMSOL model uses an exponential variogram of the form

(4)

with variance σ2 = 1 (the sill parameter) and correlation length r = 50 m (the range). For the purpose of this modeling example, assume that the parameters σ and r are known; Kitanidis (Ref. 2) gives methods for estimating their values in the case where they are unknown. Since the covariance must approach zero as the distance between samples increases, this implies a simple covariance function q(|xi−xj|) = q(h) = e−h/r.

Thus, for the penalty term, the elements of the covariance matrix Q are easy to evaluate. Computing the inverse, Q−1, would be unnecessarily expensive. Instead, split Equation 2 in two, introducing an auxiliary vector u of the same length as the unknown parameter vector s:

(5)

Solving the linear system Qu = s−Xβ for the value of u is in general cheaper than inverting Q.

In Figure 1, the eight points on the rim of the area being characterized are numbered pairwise as 1±, 2±, 3±, 4±, with plus and minus signs denoting injection wells and pumping wells, respectively.

Figure 1: Discretization of the hydraulic-conductivity parameter field. Each 10 m-by-10 m square zone is associated with a degree of freedom encoding the constant log10 Ksvalue in the zone. The numbers in blue, outside the square grid, are labels for the dipole-pump pairs (“+” for injection wells and “-” for pumping wells).

Model Data

The data required for setting up the example are supplied in text files included in your COMSOL installation:

•

A text file with reference synthetically generated field data, containing the log10 values of the regularized hydraulic-conductivity parameter field that is used to generate fictitious hydraulic-head measurements. This allows you to evaluate the optimization solver’s performance and accuracy, and to test and calibrate the inverse model. For example, you can try out different penalty terms in the objective function and investigate the solution’s dependence on the number of observations used.

The log10 Ks field for the reference model is visualized in Figure 2. It was generated with the exponential variogram, γ, given in Equation 4.

•

Four CSV-files contain the hydraulic-head measurement data series. These data sets were generated by solving the reference forward model for the four different dipole-pump configurations.

•

An auxiliary measurement data file containing a single zero value, to be used when implementing the penalty term formally on least-squares form.

Figure 2: Discretized hydraulic-conductivity field for the reference forward model.

The four simulated dipole-pump tests for the reference forward model give the combined results listed in Table 1, here rounded off to the nearest millimeter; the example uses these numbers as fictitious measurement values.

Table 1: Hydraulic Head (meters) measured at wells for the four dipole pump tests.
Test	1+	1-	2+	2-	3+	3-	4+	4-
1	-	-	1.425	-0.718	-0.028	0.333	1.197	-1.330
2	0.568	-1.576	-	-	1.425	-1.176	-0.641	-0.689
3	-0.801	-0.439	1.187	-1.414	-	-	-1.565	0.958
4	1.636	-0.891	0.501	0.453	-1.272	1.251	-	-

Although the normalization of the fitness term (see Equation 3) assumes a measurement accuracy of 1 cm, this model (in contrast to Ref. 1) adds no random error component to the above data. This is because the primary aim here is to investigate the software’s performance rather than the viability of the method for aquifer characterization.

Results and Discussion

Figure 3 and Figure 4 show the optimization results for the hydraulic conductivity obtained using only the first measurement series and all four measurement series, respectively. The corresponding values of the total mean square error are 0.318 (0.316) and 0.086 (0.087), where the values in parentheses are those obtained by Cardiff and Kitanidis (Ref. 1). Note that the error estimates are slightly different because Cardiff and Kitanidis include a synthetic measurement error of the order 1% in their observations. Figure 5 compares the results from four different inverse-modeling simulations, using 1, 2, 3, and 4 measurement series, respectively. As the figure shows, the improvement in accuracy is most marked when going from 1 to 2 pump tests (that is, from 6 to 12 observations). Including additional measurements appears to yield comparatively small benefits in accuracy. However, as discussed in Ref. 2, the posterior uncertainty in parameter estimates decreases as more measurements are added.

Figure 3: Inverse-modeling solution for the hydraulic conductivity logarithm parameter field obtained using only the first measurement series (dipole pair number 1).

Figure 4: Hydraulic conductivity obtained by inverse modeling using 24 observations.

Figure 5: Inverse-modeling results for the hydraulic conductivity field with different numbers of hydraulic-head observations taken into account: 6 (upper left), 12 (upper right), 18 (lower left), and 24 (lower right). For a color-scale legend, see Figure 2.

References

1. M. Cardiff and P.K. Kitanidis, “Efficient Solution of Nonlinear, Underdetermined Inverse Problems with a Generalized PDE Model,” Computers & Geosciences, vol. 34, pp. 1480–1491, 2008.

2. P.K. Kitanidis, “Quasi-linear Geostatistical Theory for Inversing,” Water Resources Research, vol. 31, no. 10, pp. 2411–2419, 1995.

Subsurface_Flow_Module/Fluid_Flow/aquifer_characterization

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 100.

Locate the section. In the text field, type -150.

In the text field, type -150.

Array 1 (arr1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type 3.

Locate the section. In the text field, type 100.

In the text field, type 100.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type -50,0,0,50.

In the text field, type 0,-50,50,0.

In the toolbar, click

Click the

button in the toolbar.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Description
N0	0.1[kg/(m*s)]	Pump source strength
deltaH	1[cm]	Hydraulic head measurement error
logKs0	-5	Hydraulic conductivity, initial log10 value
th	1	Measurement-series index
sigma	1	Sill parameter
r	50[m]	Range parameter
elementTypeFactor	1	1 for quads, 0.5 for triangles

As already mentioned, defining the model requires a number of external data files:

•

The hydraulic-head measurement data series are stored in four separate CSV-files included in your COMSOL installation.

•

A text file with the forward-model data for the log-transmittivity parameter field is provided to allow you to evaluate the optimization solver’s performance and accuracy, and to test and calibrate the inverse model. For example, you can try out different penalty terms in the objective function and investigate the solution’s dependence on the number of observations used.

•

A dummy measurement file containing a single zero value, in order to implement the penalty term on the same form as a least-squares error term.

To make the data for the log-transmittivity field available in the model, create an interpolation function.

logKs Reference

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file aquifer_characterization_logKs_ref.txt.

From the list, choose .

Click

In the text field, type logKs Reference.

Locate the section. Find the subsection. In the table, enter the following settings:


Function name	Position in file
logKs_ref	1

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Function	Unit
logKs_ref	1

In the table, enter the following settings:


Argument	Unit
Argument 1	m
Argument 2	m

Click

Definitions

Define Infinite Element Domains for the outer squares.

Infinite Element Domain 1 (ie1)

In the toolbar, click

Select all rectangles except the one in the middle.

Variables 1

In the window, right-click and choose .

For an initial evaluation of the forward model, use the reference solution’s hydraulic conductivity field. Later, the control variable field for the Optimization interface will override this expression.

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Description
logKs	logKs_ref(x,y)	Hydraulic conductivity, log 10 value

Materials

Assign water as the only domain material everywhere. Darcy’s law in general requires separate specifications of the fluid and matrix material properties. In this case, the matrix conductivity is given first by an interpolation field, later by the Optimization interface. Darcy’s Law for time-dependent problems also involves the porosity, but this model is stationary so you can safely ignore the warning about this property being undefined.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Darcy’s Law (dl)

Porous Matrix 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the K text field, type 10^logKs0.
This is the conductivity that is applied to the . The center domain is defined as follows:

Porous Medium 2

In the toolbar, click

and choose .

Select Domain 5 only.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the K text field, type 10^logKs.

Hydraulic Head 1

In the toolbar, click

and choose .

Select all outer boundaries.

Next, add the dipole pumps by using the feature.

Line Mass Source 1

In the toolbar, click

and choose .

Select Point 6 only.

In the window for , locate the section.

In the N0 text field, type if(th==1,N0,0).

The conditional statement allows you to control which dipole pump to activate through the parameter th. Note that a value of zero is the default condition, so setting N0 to zero has no effect on the model equations.

Line Mass Source 2-8

Proceed to create seven additional features with the following settings:


Name	Selection	Expression
Line Mass Source 2	7	if(th==2,N0,0)
Line Mass Source 3	8	if(th==3,N0,0)
Line Mass Source 4	10	if(th==4,N0,0)
Line Mass Source 5	11	if(th==4,-N0,0)
Line Mass Source 6	13	if(th==3,-N0,0)
Line Mass Source 7	14	if(th==2,-N0,0)
Line Mass Source 8	15	if(th==1,-N0,0)

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 5 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type 5.

Size 2

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Points 6–8, 10, 11, and 13–15 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type 0.1.

Click

Study 1

Solve the forward problem to verify that you have implemented the physics correctly.

In the toolbar, click

Results

Hydraulic Head, Forward

The default plot shows the pressure distribution. Modify it to display the hydraulic head instead.

In the window for , type Hydraulic Head, Forward in the text field.

Surface

In the window, expand the node, then click .

In the window for , click in the upper-right corner of the section. From the menu, choose .

The result should look like that in the figure below.

Next, set up the auxiliary model for the hydraulic conductivity control variable field and the help variable u.

Add Component

In the window, right-click the root node and choose .

Next, create the auxiliary geometry on which the Optimization interface lives. This geometry is a copy of the aquifer domain, except that you do not need to include vertices at the positions of the sources and sinks.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the window, click the root node.

In the toolbar, click

to close the window.

Geometry 2

In the window, under click .

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 100.

Locate the section. From the list, choose .

Definitions (comp2)

You will define the hydraulic conductivity control variable shortly. To make it available on the forward-model’s geometry, use a General Extrusion operator.

General Extrusion 1 (genext1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Define an operator for calculating the mean square error. Note that this is an area-weighted average.

Average 1 (aveop1)

In the toolbar, click

and choose .

In the window for , type mean in the text field.

Locate the section. From the list, choose .

To evaluate the penalty function, which is defined in terms of matrix-vector products with discrete control variable degrees of freedom, you need a way to compute a discrete sum over the mesh elements. This can be achieved using an Integration operator of order 0, together with appropriate compensation for the mesh element area.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , type int0 in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type 0.

Variables, Domain

In the toolbar, click

In the window for , type Variables, Domain in the text field.

Locate the section. From the list, choose .

Select Domain 1 only.

Locate the section. In the table, enter the following settings:


Name	Expression	Description
logKs_ref	logKs_ref(x,y)	Hydraulic conductivity, reference model
areaFactor	1/(elementTypeFactor*dvol)	Summation compensation factor
dist	sqrt((x-dest(x))^2+(y-dest(y))^2)	Distance between points inside summation

The dest() operator instructs the software to evaluate the enclosed argument at the destination point when evaluating an integral. The result of calling a nonlocal integration coupling with an argument containing the variable dist will therefore depend on where the integral is evaluated.

In this case, a nonlocal integration coupling is used to compute a discrete sum over elements. This is achieved using a single integration point per element and compensating for the corresponding integration point weight. The variable dvol is the volume scale factor that the software uses internally when mapping between dimensionless local element coordinates and the model geometry’s global coordinate system shown in the user interface. The mesh-element volume in the local element coordinate system multiplied by dvol therefore gives the mesh-element volume in the global coordinates. In 2D, the mesh-element area in element coordinates is 1 for quadrilateral meshes and 1/2 for triangular meshes.

Variables, Global

In the toolbar, click

Add variables for the mean squared error, the mean of the log conductivity and the penalty term.

In the window for , type Variables, Global in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Description
MSE	mean((logKs-logKs_ref)^2)	Area-weighted mean squared error
logKs_mean	int0(logKs*areaFactor)/int0(areaFactor)	Discrete control variable mean
L_penalty	int0((logKs-logKs_mean)uareaFactor)	Penalty function

Note that the penalty function is evaluated using the intermediate variable u, which is computed as the solution of an auxiliary equation. The integration operator int0 together with areaFactor effectively computes the discrete scalar product between the vectors logKs-logKs_mean and u.
Furthermore, note that the model entry appears in red as logKs has not yet been defined in component 2. This will be done in the General Optimization settings.

Covariance function

In the toolbar, click

In the window for , type Covariance function in the text field.

In the text field, type Q.

Locate the section. In the text field, type sigma^2*exp(-x/r).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	m

In the text field, type 1.

General Optimization (opt)

In the window, under click .

Control Variable Field 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type logKs.

In the text field, type logKs_ref.

This is just a temporary initial setting; after creating a plot of the reference solution you will change it to logKs0.

Locate the section. From the list, choose .

Find the subsection. From the list, choose .

This gives discontinuous elements of order zero for the control variable logKs, a choice that ensures that each auxiliary-geometry mesh element carries a single degree of freedom.

Global Least-Squares Objective 1

Here, use the measurement file aquifer_characterization_zero.csv to minimize a positive objective contribution that does not correspond to an actual measurement; you will add the real measurements shortly to the first model.

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file aquifer_characterization_zero.csv.

To put the penalty term formally on least-squares format — something which is required by the Levenberg-Marquardt solver — enter the square root of the penalty value as a model expression to be compared with the dummy zero measurement.

Value Column 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type sqrt(L_penalty).

Domain ODEs and DAEs (dode)

Set up the equation for the auxiliary vector u, which must use the same discretization as the control variable field logKs.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Source term quantity	Unit
Custom unit	1

Click to expand the section. From the list, choose .

Distributed ODE 1

In the window, under click .

In the window for , locate the section.

In the f text field, type (logKs-logKs_mean-int0(u*Q(dist)*areaFactor)).

Mesh 2

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 10.

Click

Definitions (comp1)

Now return to and add an alternative definition of variable logKs, mapping the corresponding control variable field from .

Variables 4

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Description
logKs	comp2.genext1(comp2.logKs)	Hydraulic conductivity, log 10 value

Component 1 (comp1)

Because the measurements are to be compared with the solutions of , define the least-squares objective contributions in this model.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

General Optimization 2 (opt2)

Least-Squares Objective 1

Right-click and choose the domain setting .

Select Domain 5 only.

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file aquifer_characterization_H1.csv.

This file has three columns: the first two are coordinate columns defining the measurement location for the value found in the third column. The index 1 in the filename refers to the measurement series number. Thus, this feature contains the contribution from the first measurement series to the objective function, corresponding to the experimental parameter value th=1. Once you have set it up, you can duplicate it to conveniently add the contributions from the remaining three experiment series.

Locate the section. Click

In the table, enter the following settings:


Name	Expression
th	1

Coordinate Column 1

In the toolbar, click

and choose .

Least-Squares Objective 1

In the window, click .

Coordinate Column 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Least-Squares Objective 1

In the window, click .

Value Column 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type comp1.dl.H.

In the text field, type 1/deltaH^2.

The measurements are weighted by the inverse of the square of the measurement error.

Least-Squares Objectives 2-4

Now add the three remaining measurement series, by duplicating the node for each measurement. Make sure you end up with the following nodes and the correct settings:


Name	Filename	Parameter th
Least-Squares Objective 2	aquifer_characterization_H2.csv	2
Least-Squares Objective 3	aquifer_characterization_H3.csv	3
Least-Squares Objective 4	aquifer_characterization_H4.csv	4