Adjoint Method for an Inverse Problem of CCPF Model∗

Zhenhua CHENKaiqi ANYuan LIUWenbin CHEN

1 Introduction

Groundwater,being one of the most useful water sources in worldwide human society,meets one-third of the overall water needs in France and supplies up to 10%–15%of the water consumption in China.Among several different types of groundwater,karst aquifer should be mentioned here as one typical system.It mainly consists of the porous medium,referring to the matrix,which serves to hold the water.And in this system,some water conduits or pipes are embedded among the matrix which plays an important role in transporting the fluid flow and the contaminates.

For modeling the karst aquifers,the most used approach for geological studies is referred to as the coupled continuum pipe flow(CCPF,for short)model.The model is a coupled system consisting of a two-(2D,for short)or three-dimensional(3D,for short)continuum domain Ωm,namely the matrix.The flow in matrix is governed by a Darcy-type system.Inside the matrix lies a one-dimensional(1D,for short)conduit Ωc,which is governed by the pipe flow model.It was used in[2,4,15]to study the gensis of karst aquifers firstly.And the model,to some extent,reaches success.In the literature[14],a modified continuous CCPF model was proposed,which was proved to be well-posed in the 2D case.Furthermore,Wang studied the well-posedness of the model in 3D case(see[19]).

In the CCPF model,the fluid exchange between the conduit and the matrix is calculated via an ad-hoc termα(hm−hc),whereαis the exchange rate,andhmandhcare the hydraulic heads in the matrix and conduit,respectively.Chen et al pointed out that as long as the exchange rate was well selected,the relative error between the CCPF model and the Stokes-Darcy system could be less than 1%(see[9]),which implied the validity of the CCPF model in describing the fluid flows in karst aquifers as well as the significance of the exchange rateα.Since the discharge of the exchange flow varies in different locations,a precise description of the model entailsαto be a function that depends on the horizontal variablex.For the forward problem,the well-posedness and regularity of the weak solution have been obtained in[8,14,17,19].And for the inverse problem,in[16],Lu et al firstly verified the uniqueness of the exchange rate and proposed an inverse problem of determining the exchange rate functionα.

Usually,the inverse problem is often turned into one optimization problem and solved by the gradient-based method,which is generally considered to be accurate but inefficient(see[6,18]).In this paper,we regard the inverse problem as a PDE-constraint optimization problem where the exchange rate function minimizes a functional that measures the error between the computed and measured data(see[13]).The solution of this constrained optimization problem is then accomplished by one gradient-based optimization algorithm that requires the computation of the derivative with satisfaction of the constraints.An auxiliary problem,called the adjoint problem,is solved in order to calculate the derivative efficiently.The so-called adjoint method has successfully applied in different fields(see[7,11–12]).The forward and adjoint problems are solved by using efficient finite difference methods and iterative solvers.Similar technique has been used to reconstruct the interface of discontinuity in the conductivity(see[10]).

At each iteration of the optimization algorithm,the main computational cost of the derivative calculation is independent of the number of parameters used to represent theαsince only the solutions of a forward problem and an adjoint problem are needed,a new search direction is calculated by the steepest descent update,and the exchange rate function is updated at each iteration by new search direction.

The rest of this paper is organized as follows.In Section 2,the framework of adjoint method for PDE-constraint optimization problem is presented.In Section 3,the froward problem and inverse problem of the CCPF model are introduced,and the adjoint gradient representation calculation and algorithm for the inverse problem are given in Section 4.Finally,numerical results verify the algorithm in Section 5.

2 Framework of Adjoint Method for PDE-Constraint Optimization Problem

In this section,the framework that will be used for reduced derivative computation is presented in a functional analytical setting(see[13]).We first transform the general optimization problem into a reduced optimization problem,then the optimality conditions and an adjoint based representation for the reduced gradient of the objective function are stated.

2.1 Optimization problem

Formally,we will consider one problem having the following form:

whereJ:W→R is the objective function,E:W→Zandc:W→Rare operators,W,Z,Rare real Banach spaces,K⊂Ris a closed convex cone,C⊂Wis a closed convex set.

In general,W,ZandRare generalized function spaces and the operator equationE(w)=0 stands for a PDE or a system of coupled PDEs,constraintc(w)∈Kis regarded as an abstract inequality constraint.In some cases,it is convenient to writec(w)∈Kas the form ofw∈C,whereC⊂Wis a closed convex set,as a result,the inequality constraint is dropped:

For PDE-constraint problem,we will have one additional structure:Optimization variablewcould be divided into two parts,a statey∈Yand a control(or design)u∈U,whereYandUare Banach spaces.ThenW=Y×U,w=(y,u),and the problem becomes

Herey∈Yis described byE(y,u)=0(usually a PDE).The control(or design)u∈Uis a parameter that shall be adapted in an optimal way.

2.2 Reduced problem and adjoint method

We consider again the optimal control(or design)problem of the form

whereJ:Y×U→R is the objective function,E:Y×U→Zis an operator between Banach spaces,Wad⊂W:=Y×Uis a nonempty closed set.Assuming thatJandEare continuously Fr´echet-differentiable,and for anyu∈U,the state equationE(y,u)=0 has unique solutiony(u)∈Y.Thus,we have a solution operatoru∈U 7→y(u)∈Y.Furthermore,we assume thatEy(y(u),u)∈L(Y,Z)is continuously invertible,then the implicit function theorem ensures thaty(u)is continuously differentiable.By differentiatingE(y(u),u)=0 with respect tou,we obtain an equation fory′(u):

After insertingy(u)into(2.4),we get one reduced problem:

For the optimal problem,the computation of the derivative of the reduced objective functionis very important.There are two methods to do this:one is the sensitivity approach,the other is the adjoint approach.For the sensitivity approach,the computational cost grows linearly with the dimension ofU(see[13]).

2.2.1 Adjoint approach

In fact,there is one efficient way to compute the derivative of.From

we see thatConsequently,the requirement of the operatory′(u)∈L(U,Y)can be weakened,and only the vectory′(u)∗Jy(y(u),u)∈U∗is really required.By(2.5),y′(u)∗Jy(y(u),u)=−Eu(y(u),u)∗Ey(y(u),u)−∗Jy(y(u),u),and by introducing an adjoint stateλ=λ(u)∈Z∗,whereλ(u)solves the adjoint equation

theny′(u)∗Jy(y(u),u)=Eu(y(u),u)∗λ(u).So we haveThe derivativecould be computed by the adjoint approach as follows:

1.For givenu,we solve the state equationE(y(u),u)=0,then we get the state variabley(u).

2.By solving the adjoint equation(2.7),we get the adjoint variableλ(u).

3.Compute

which is the adjoint gradient representation.

2.2.2 A Lagrangian-based view of the adjoint approach

The adjoint gradient representation could be derived in a more general way.For optimization problem(2.4),let us define Lagrange functionL:Y×U×Z∗→R,

Substitutingy=y(u)gives,for arbitaryλ∈Z∗,

And differentiating the above equation gives

One specialλ=λ(u)can be chosen such that

This is the adjoint equation,actually,for one variational directiond∈Y,

thenLy(y(u),u,λ)=Jy(y(u),u)+Ey(y(u),u)∗λ.According to(2.9),we obtain that

By solvingE(y,u)=0,we gety.Furthermore by solving adjoint equationEy(y(u),u)∗λ=−Jy(y(u),u),we getλ.By(2.8),the adjoint gradient representation reads

2.3 Optimality conditions

Let us consider the problem(2.4).By the definition ofUadin(2.6),the optimization problem is equal toJ(y,u)s.t.E(y,u)=0,u∈Uad.We introduce the following assumption.

Assumption 2.1Assume that

1.Uad⊂Uis nonempty,convex and closed.

2.J:Y×U→R andE:Y×U→Zare continuously Fr´echet-differentiable,andU,Y,Zare Banach spaces.

3.For allthere exists a neighborhoodandE(y,u)=0 has a unique solutiony=y(u)∈Y.

4.For allderivativeEy(y(u),u)∈L(Y,Z)is continuously invertible.

For the reduced problem(2.6),we have the following general result(see[7,13]).

Theorem 2.1Assume Assumption2.1holds.Ifis a local solution of the reduced problem(2.6),thensatis fies the variational inequality

Next,we use the adjoint representation of derivative(2.10)

whereλ(u)∈Z∗is the solution of the adjoint equationEy(y(u),u)∗λ=−Jy(y(u),u).Recalling the Lagrange function associated with(2.6)L:Y×U×Z∗→R,

we get the corollary of Theorem 2.1.

Corollary 2.1Let(,)be an optimal solution of(2.6)and assume Assumption2.1holds.Then there exists an adjoint state(or Lagrange multiplier)∈Z∗,such that the following conditions hold:

Using Lagrange function,the above equations could be written in the compact form:

In this paper,we are interested in the numerical efficiency of the adjoint approach to solve the inverse problem,and do not want to verify Assumption 2.1,which can be verified in[13].

3 Analysis on the Mathematical Model

3.1 CCPF model

CCPF model is a coupled system consisting of a 2D or 3D continuum domain Ωm,which stands for the matrix.There exists a 1D conduit Ωcin the matrix,and the situation in the conduit is governed by pipe flow model.In practice,the matrix is soil or rock aquifers,the fluid exchanges between the matrix and the conduit(see Figure 1).In this paper,we assume that the matrix domain is Ωm:=(0,L)×(−M,M),Ωcis the conduit domain,andνdenotes the tangential direction along the 1D pipe/channel conduit.

Figure 1 The geological description of the karst aquifers.hmis the hydraulic head in the matrix,and hcis the hydraulic head in the conduit

We first review the formulation of the coupled pipe flow/Darcy model.The flow in the porous matrix is modeled by Boussinesq equation(see[3]):

wherehmdenotes the hydraulic head in the porous matrix,K∈R2×2is a hydraulic conductivity tensor.Sis the storativity,andγrepresents the volumetric rate offluid transfering to the porous matrix from the conduit system per unit length.

For conduit flow,the discharge is related to the head difference in the tube by applying the Darcy-Weisbach equation(see[5])

wherehcis the hydraulic head in the conduit,Poisieulle constantheredis the pipe diameter(channel width for 2D),Qis the total discharge in the pipe.Conversation of mass in the pipe impliesthen we haveThe matrix and conduit flows are coupled at their intersection by the quasi-steady-state exchange term(see[8,17,19]),

whereα>0 is the exchange rate function,this means that the process in the conduit is enslaved by that in the porous matrix.

As a result,in the steady-state case,we have the system

wherehmis the hydraulic head in matrix,hcis the hydraulic head in conduit,δΩcis Diracδfunction concentrated on Ωc,fmandfcrepresent the external source or sink.Laminar Poisieulle constantwheredis the diameter of the conduit,gis the earth’s gravitational acceleration,andρis the kinematic viscosity of water,is the exchange rate function wheresdenotes arc length variable.

3.2 Well-Posedness and regularity

Under the case of Dirichlet boundary con dition, fixed hydraulic heads on the boundaries in both domains are introduced by

wheregmandgcare given a priori.

Define the spaceWe introduce a bilinear form as follows:

whereh=(hm,hc)∈H,and the test functionsv=(hm,hc)∈H.Then for allh∈H,the weak solution of(3.4)hyields

The following theorem shows that the weak solution uniquely exists under appropriate assumptions.

Theorem 3.1Assume thatThen theweak solution of(3.4)huniquely exists,satisfying the following estimation:

where C is a constant independent offmand fc.

Higher regularity of the solution is possible,the proof of the theorem can be found in[8,17,19].The well-posedness of the weak solution for Neumann boundary conditions is similarly compared with the Dirichlet boundary conditions.

3.3 Inverse problem of determining the exchange rate function α

In CCPF model,the fluid exchange rate between the matrix and the conduit is calculated by an ad-hoc termα(hm|Ωc−hc),whereα:=α(s)is the exchange function de fined on the conduit domain Ωc,andsstands for arc parameter.For anisotropic exchange rate functionTheorem 3.1 ensures the well-posedness of CCPF model.Now we discuss the uniqueness of the exchange rate by measuring the Cauchy data along one side of the boundary.

In the numerical realization,the geometry is described as following:The matrix Ωmis a rectangle domain,Ωm:=(0,L)×(−M,M),and the conduit Ωccould be expressed byy=ψ(x),whereψ(x)∈C2([0,L])andψ(0)=ψ(L)=0.Furthermore,we assume that the choice ofψ(x)allows the well-posedness of the elliptic problem in a Lipschitz subdomain,wheresatis fiesLet Γ:={0}×(−M,M).The Cauchy data of partial boundaries of conduit domain Ωcand matrix domian Ωmare measured:

Theorem 3.2Let us define k(x)=α(s(x))(hm(x,ψ(x))−hc(s(x))).If there exist two functions k1(x)and k2(x)satisfying CCPF model

with the same Cauchy data(3.8)–(3.11),there holds k1(x)=k2(x)almost everywhere.

The proof of the theorem can be found in[16].In addition,ifk(x)=α(s(x))(hm(x,ψ(x))−hc(s(x))),let us define

Under some assumptions,if there exist two exchange rate functionsα1(s)andα2(s)having the same Cauchy data(3.8)–(3.11),there holdsα1(s)=α2(s)almost everywhere(see[16]).

Then the inverse problem could be expressed as follows:if the Cauchy data are given on partial boundaries of matrix and conduit,we try to find an exchange rate functionα(s)such that the associated weak solutionhequals to the Cauchy data on the partial boundary.

4 Adjoint Method for Inverse Problem of CCPF Model

The inverse problem of CCPF model can be posed as optimization problem with the constraint given by a variational equation.The cost functional is the mismatch between the state obtained from a trial solution and the measured one:

where Ω1and Ω2are the parts of∂Ωmand∂Ωc,respectively.Letbe the measured data andh=(hm,hc)the weak solution of(3.4).For instance,let us consider the Neumann boundary conditions:

where Γ1:=L×(−M,M)is the right boundary of Ωm,andSis the length of Ωc(see Figure 2).

Figure 2 Domain and boundaries.

Usually,we consider the no flow boundary condition aty=±M,

The physical background for this particular Neumann condition is that we assume there is no flows through the matrix on the upper and lower boundaries.

De fine the space:=H1(Ωm)×H1(Ωc).The weak formulation is now: findingh∈such that

which is one PDE-constraint.Consequently,solving the inverse problem equals to solving an optimization problem with the PDE-constraint,

Define the spaces

Then the weak formulation(4.5)is equivalent to: findingh∈such that

The formulation now de fines the state equation operator

Thus,we can de fine the reduced problem

4.1 Calculation of the gradient

From the adjoint gradient representation(2.10),we have to compute the Lagrange multiplierλ=(λm,λc)∈Z∗by solving the adjoint system(2.7),which reads as

with the given weak solutionhof the state equation.The partial derivativeEhwith respect tohis Fr´echet derivative with the direction ofω=(ωm,ωc).For the forward problem(4.8),the adjoint equation is:for allω∈,

where

If∆αis the variational direction of the derivative,then

We consider an objective functional of the type

If∆αis the variaitonal direction ofα,then

Finally,we get the adjoint gradient representation

4.2 Optimization algorithm

Practically,the optimal variableα(s)is defined by the design parameteru∈Uwith a finite or infinite dimensional design spaceU.Thus,we have a mapα:u→α(u).Using the chain rule,derivatives of the objective functionare obtained as

If the design space ofα(s)is finite dimensional,where(ϕ1,···,ϕd)is a basis of design spaceU,and thenwithwherehandλare solutions of the state equation and the adjoint equation,respectively,Gi=G(h,λ ,ϕi)is the element of the gradient,and will be supplied to an optimization algorithm which is used in the solution of the inverse problem.For example,for the steepest decent method,we have the following algorithm.

Algorithm 4.1(Steepest Decent Method)Choose an initial vectorassociated with the initial guessα0(s).

Fork=0,1,2,···,

1.IfSTOP.

2.Choose a descent directiondk∈U,such that

3.Choose a step sizeσk>0 such that

4.Set

In our numerical experiments,the step size rule we used here is the Armijo rule.

Armijo RuleGiven a descent direction∆αof,choose the maximumfor which

whereγ∈(0,1)is a constant.

On the other hand,the descent direction∆αcould be chosen by a special rule.

Lemma 4.1If A is a Hilbert space,b(·,·)is a uniformly bounded and coercive bilinear form in A that satisfies

where α0,β0>0,and∆α satisfies

then it is guaranteed that ζ∆α decreases the cost functional J,where ζ is a small enough positive number.

ProofFrom the definition,for any∆α∈A,

and

Indeed,the second term on the right side ofinequality is strictly negative and the third term can be chosen such that

Theorem 4.1(Choice of the Descent Direction)Ifh=(hm,hc)is the weak solution of(4.8),λ=(λm,λc)is the weak solution of adjoint equation(4.10),then we can choose

such that

where ζ is small enough.Choice of ζ could use the Armijo rule.

ProofFrom Lemma 4.1,the bilinear form isand obviously it satisfies(4.16)–(4.17).The adjoint gradient representation gives

Because of the choice of the∆α,

we get

5 Numerical Results

In this section,we demonstrate the adjoint derivative calculus on a numerical model problem.We choose the decent direction∆α(s)by Theorem 4.1.Numerical results show the efficiency.

5.1 Problem description

For convenience and simplicity,the computational domain in the model is the matrix domainΩm:=[0,1]×[−1,1]with the upper sectionand the lower section[0,1]×[−1,0],and the conduit domain is Ωc=[0,1]×{0}(see Figure 3).

We begin to construct one exact solution of the model with a target exchange rate function,for convenience of numerical experiments,the hydraulic conductivity tensor K is set to be I,and the Poisieulle constantDis set to be 1.Furthermore,considering the fact that the domainΩcis a straight line in the experiments,the CCPF model in this section is

We adjust the forcingfmandfcin Ωmand Ωc,respectively,such that the exact solution is

Figure 3 Computational domain in the numerical examples.

As a result,the target exchange rate function is 2+sin(πx),fmandfcin(5.1)are

5.2 Convergence of the finite difference method

For the forward problem and the adjoint problem,we use finite difference method to solve them.Here we present the convergence of the finite difference method,i.e.,we will prove that the numerical approach is accurate for solving the partial differential equation systems.We compute the finite difference approximation by using the sequence of grid sizesh=2k,k=1,···,6.Errors are measured by the discretel2norm.

Definition 5.1The convergence factor is given as convergence factorwherehkdenotes the solution under grid size2k.

In Figure 4,the convergence factor is plotted versusk.From the figure,we know that the solver to the forward problem has second order accuracy.

5.3 Results ofinverse problem

If the regular termis ignored,we recast the problem as a PDE-constraint optimization problem

whereandare the measured data,i.e.,the exact solution of(5.1)with the target exchange rate function 2+sin(πx).At every step of the iteration,∆αis chosen by Theorem 4.1.And the initial guess ofαisy=2.The iteration stops when

Figure 4 Convergence factor offinite difference approximation of h versus k for forward problem.

or the number of the iteration reaches the maximal step,whereJ(i)denotes the energy of stepi.The numerical results are sorted according to Ω1and Ω2.

Case 1:Ω1={0}×[−1,1],Ω2={0},i.e.,measured data on the left side of the boundary are given.Boundary conditions:

Case 2:Ω1={0}×[−1,1]∪{1}×[−1,1],Ω2={0}∪{1},i.e.,measured data on the left and right sides of the boundary are given.Boundary conditions:

Case 3:For ideal case,Ω1=Ωm,Ω2=Ωc,i.e.,all the measured data are given and the boundary condition is Dirichlet boundary condition.

Figures 5–7 correspond to Cases 1–3,respectively.Comparing these three results,we can find more accuracy of the unknown value if we know more information of the solutions.Specially for the Case 3 where the solutions are known in whole domain.Based on the results from Figure 7,it is hopeful to get one theoretical convergence analysis for the unknownα.But if only partially data(such as boundary values of the solutions)are measured as Case 1 and Case 2,the problems turn to be ill-posed,and it is hard to get enough accuracy.

Next,we discuss the case that the measured data are perturbed by noise.For the above three cases,we add±5%noise to the measure data,the numerical results are reported in Figure 8(Case 1),Figure 9(Case 2)and Figure 10(Case 3).In the cases with noises,these figures show that the energy functional stagnates in the higher level than the cases without noises,while the inversion results can still be acceptable.

Figure 5 Iteration process(left)and log of the energy(right),k=5,measured data of left side are known,iteration number is 20700.

Figure 6 Iteration process(left)and log of the energy(right),k=5,measured data of left and right sides are known,iteration number is 592.

Figure 7 Iteration process(left)and log of the energy(right),k=5,measured data of all the domain are known,iteration number is 40000.

Figure 8 k=4,measured data of left side are known,noise level:±5%,iteration number is 1098.

Figure 9 k=4,measured data of left and right sides are known,noise level:±5%,iteration number is 1098.

Figure 10 k=4,measured data of all the domain are known,noise level:±5%,iteration number is 550.

AcknowledgmentsThe authors thank Xiaoming Wang,Neil Hua and Shuai Lu for their helpful discussions and suggestions.

[1]Adams,R.A.,Sobolev Spaces,Academic Press,New York,1975.

[2]Bauer,S.,Liedl,R.and Sauter,M.,Modeling of karst aquifer genesis:Influence of exchange flow,Water Resour.Res.,39(10),2003,SBH6.1–SBH6.12.

[3]Bear,J.and Verrruijt,A.,Modeling Groundwater Flow and Pollution,D.Reidel.Pub.Co.,Norwell,Mass,1987.

[4]Birk,S.,Liedl,R.,Sauter,M.and Teutsch,G.,Hydraulic boundary conditions as a controllong factor in karst gensis:A numerical modeling study on artesian conduit development in gypsum,Water Resour.Res.,39(1),2003,SBH2.1–SBH2.14.

[5]Bobok,E.,Fluid Mechanics for Petroleum Engineers,Elservier Sci,New York,1993.

[6]Bonnans,J.F.,Gilbert,J.C.,Lemarechal,C.and Sagastizabal,C.A.,Numerical Optimization:Theoretical and Practical Aspects,2nd Edition,Springer-Verlag,Berlin,2006.

[7]Brandenburg,C.,Lindemann,F.,Ulbrich,M.and Ulbrich,S.,A continuous adjoint approach to shape optimization for Navier-Stokes flow,Optimal Control of Coupled Systems of Partial Differential Equations,International Series of Numerical Mathematics,158,2009,35–56.

[8]Cao,Y.,Gunzburger,M.,Hua F.and Wang X.,Analysis and finite element approximation of a coupled,continuum pipe-flow/Darcy model for flow in porous media with embedded conduits,Numer.Math.PDE,27(5),2010,1243–1252.

[9]Chen,N.,Gunzburger,M.,Hu,B.,et al.,Calibrating the exchange coefficient in the modified coupled continuum pipe-flow model for flows in karst aquifers,J.Hydrol,414–415,2012,294–301.

[10]Chen,W.,Cheng,J.,Lin,J.and Wang,L.,A level set method to reconstruct the interface of discontinuity in the conductivity,Science in China Series A:Mathematics,52(1),2009,29–44.

[11]Feij´oo,G.R.,Malhotra,M.,Oberai,A.A.,et al.,Shape sensitivity calculations for exterior acoustics problems,Eng.Comput.,18(3/4),2001,376–391.

[12]Feij´oo,G.R.,Oberai,A.A.and Pinsky,P.M.,An application of shape optimization in the solution ofinverse acoustic scattering problems,Inverse Problems,20(1),2004,199–228.

[13]Hinze,M.,Pinnau,R.,Ulbrich,M.and Ulbrich,S.,Optimization with PDE Constraints,Springer-Verlag,Berlin,2009.

[14]Hua,F.,Modeling,analysis and simulation of Stokes-Darcy system with Beavers-Joseph interface condition,Ph.D.Thesis,Florida State University,Tallahassee,2009.

[15]Liedl,R.,Sauter,M.,Hckinghaus,D.,et al.,Simulation of the development of Karst aquifers using a coupled continuum pipe flow model,Water Resour.Res.,39(3),2003,SBH6.1–SBH6.11.

[16]Lu,S.,Chen,N.,Hu,B.and Cheng,J.,On the inverse problems for the coupled continuum pipe flow model for flows in karst aquifers,Inverse Problems,28(6),2012,065003.

[17]Narasimhan,T.N.,Multidimensional numerical simulation offluid flow in fractured porous media,Water Resour.Res.,18(4),1982,1235–1247.

[18]Vogel,C.R.,Computational methods for inverse problems,Frontiers in Applied Mathematics,SIAM,Philadelphia,2002.

[19]Wang,X.,On the coupled continuum pipe flow model(CCPF)for flows in Karst aquifer,DCDS-B,13(2),2010,489–501.