Xinming WUPhilipp KÜGLERShuai LU

1 Introduction

Modeling of karst aquifers has attracted great interest nowadays because ofits important role as groundwater resources that are increasingly contaminated by industrial accidents and human settlements.Geologically,this groundwater system shall consist of a fissured volume(matrix)with high storage and low hydraulic conductivity,i.e.,limestones or dolomites,and karst channels(conduits)characterized by fast transport of water.Both structures,revealed by the geologists,own two different types offlow,for instance the diffusive flow in the matrix and turbulent flow in the conduit.Another important ingredient in describing these complex nested structures is the exchange effort between this dual flow system controlled by the difference in hydraulic heads(see Figure 1).For further details we refer to[16]and references therein.These observations promote varieties of models,for instance the Navier-Stokes/Darcy system(see[8]and references therein),the Stokes-Brinkman system(see a recent publication[7])and the coupled continuum pipe flow model[2–4,13,18,21].

Figure 1 Exchange effort between the matrix and the conduit controlled by the difference in hydraulic head.hmrepresents the hydraulic head in the matrix and hcin the conduit,respectively(see[11]).

The focus of the current work is on the CCPF model which is ad-hoc arising from the karst aquifer genesis(see[2–4]).In this particular model(2D setting)conduits are degenerated into 1D traces embedded in a 2D matrix and divided into segments with a set offinite nodes{xi}(see Figure 2).The Darcian flow is treated as a continuous flow field in the matrix Ωmby

where hmis the hydraulic head,K is the hydraulic conductivity,S is the storativity coefficient in the non-steady case,fmis the recharge rate in the continuum matrix and Πexis the most important term characterizing the exchange effort between the matrix and the conduit.The mathematical formulation of the exchange rate Πexis given by

and

where δ is the Dirac delta function,V is the unit volume of the continuum matrix and qex,iis the Barenblatt type exchange flow at the node xiwhich is assumed to be proportional to the hydraulic head difference between the matrix(i.e.,hm,i)and the conduit(i.e.,hc,i).The so-called exchange coefficient αex,iis crucial and depends on many aspects,for instance,the exchange surface,the hydraulic conductivity in the matrix and the local conduit geometry(see[2]).At the same time,by assuming that the conduit flow is laminar and obeys the Kirchhoff’s rule,the governing equation for the head hcin the conduit Ωcis concentrated at the nodes{xi}of the segments by

where fc,istands for the recharge rate to the conduit at the ith node,and the Poiseuille flow formula Qijtakes the form

with a laminar Poiseuille constant D and a segment length Lijbetween the ith and jth nodes.

Figure 2 Conduits are degenerated into 1D traces embedded in a 2D matrix and divided into segments with nodes{xi}.

The original CCPF model(1.1)–(1.2)is intuitively natural and clear but mathematically ill-posed.[21]showed that the hydraulic head hmblows up at the nodeby a Green function representation.To overcome the point singularity in the original CCPF model a series of modifications was carried out in[5,13,21]in 2D and 3D setting,respectively.In the considered scenario,the following modified steady CCPF model is proposed in the dissertation(see[13]):

where α is the exchange rate function depending on the space variable x,denotes the tangential derivative along the conduit and δΩcis the Dirac delta function concentrated on Ωc.Well-posedness of the modified version was provided in[13],its FEM numerical approximation was given in[5].The heuristic derivation from the original CCPF model to the modified one was well addressed in[21]where the exchange rate function α in(1.3)is adjusted accordingly.In principle the exchange coefficient αexin the original CCPF model(1.1)–(1.2)must differ from the one in the modified model(1.3)depending on different conduit variables.

Calibration of the exchange coefficient αexor rate function α in both CCPF models is of great importance as it determines exchange effects between the fissured matrix and the karst conduits.

In early literatures[16]provided comprehensive numerical illustration by tuning various values of the source term in(1.1)as exchange effort with different conduit discharge,i.e.,natural springs or drainage galleries.[2–3]investigated karst aquifer genesis by the original CCPF model with different constant values of αexin order to obtain different geological developments.Recently[6]focused on the validity of the modified CCPF model by calibrating a constant exchange parameter α in order to fit the Stokes-Darcy system which is viewed as a“true model”.Nevertheless,calibrating the exchange rate function directly is difficult in real situations because the conduit is embedded in the matrix and usually not easy to reach.Thus,to identify this function from indirect measurements is ofinterests.However,to reconstruct such a function defined on an interior trace in the CCPF model is rather new and,to the best of our knowledge,only[17]reported a uniqueness result from Cauchy data on parts of the boundary,for instance,by measuring the hydrology head and the seepage velocity.No algorithms have been proposed to numerically approximate α from indirect measurements.

At the same time,the reconstruction of the exchange rate function α from indirect data can be viewed as a class of nonlinear inverse problems that has been well-investigated from different aspects including the geophysical identification or industrial applications(see[9]).Often,such nonlinear inverse problems are ill-posed,for which compactness of the forward mapping together with,e.g.,its(local)injectivity is a sufficient criterion(see[9]).In order to solve ill-posed problems,we require appropriate regularization techniques.In Section 3 a short introduction with two iterative regularization schemes will be provided as background knowledge.A recent monograph on the subject is[15]which contains comprehensive convergence analysis for these approaches.We note that in order to numerically demonstrate the procedure,we will confine our problem to a rectangle domain and a straight conduit,but a generalized setting can be established in a similar manner.

The paper is organized in the following sense.In Section 2 we formulate an abstract parameter-to-output nonlinear operator which maps the exchange rate function α to the local Neumann boundary data(i.e.,the measurable seepage velocity of the fluid flow).The mathematical compactness of this operator is verified as well.The corresponding inverse problem ofidentifying the exchange rate function from the measured(noisy)boundary data and its iterative regularization schemes are thus considered in Section 3 in order to obtain a stable approximation.Section 4 contains several numerical examples which verifies the appropriateness of our proposed approaches.

2 Well-Posedness of the Direct Problem

As mentioned in the previous section,the CCPF model in the current work will be considered in a simplified geometry of the domain with homogeneous,the isotropic porous media and a straight conduit.To be more precise,the matrix and conduit are defined as Ωm=(0,1)×(−1,1)and Ωc=(0,1)×{y=0},respectively;the hydraulic conductivity K is equal to KI with a constant K and an identity tensor I;D represents a Poiseuille constant.We thus obtain the following CCPF model for a laminar flow:

with Dirichlet boundary conditions

Assumption 2.1In the sequel of this paper,we will assume that the source terms satisfyfm∈L2(Ωm)andfc∈L2(Ωc).Moreover,the exchange rate functionαbelongs to a closed subsetL2(Ωc)∩Π whose definition is

Remark 2.1The well-posedness of the CCPF model(2.1)–(2.2),i.e.,existence,stability as well as uniqueness of the weak solution,can be proven by assuming weaker regularitiesα(Ωc),fm∈H−1(Ωm)andfc∈H−1(Ωc)referring to[5,13].Nevertheless,higher regularities listed in Assumption 2.1 are necessary to well define a so-called parameter-to-output mapping(operator)Fin the forthcoming equation(2.6).

Under Assumption 2.1 the following estimates hold with homogeneous Dirichlet boundary conditions respectively for(2.1):

with different regularity of the right-hand sides.We note that the global regularity ofhmin(2.4)is nearly optimal in the sense that there exists no constantβ≥0 satisfyinghmeven assuming thatfmandfcare smooth.For detailed arguments we refer to[5,13].

In order to put the following contents into the classical framework of regularization theory ofinverse problems(see[9]),we will formulate the inverse problem as an abstract(nonlinear)operator equation,i.e.,

whereαis the exchange rate function andzrepresents the observed data.Referring to Section 1,the observed data can be chosen as the seepage velocity of the fluid flowK∇hmat the boundaries∂Ωm,∂Ωcwhich is equivalent to partial Neumann data of(hm,hc).Without loss of generality,by denoting Γ1={x=0,y∈(−1,0)},Γ2={x=0,y∈(0,1)}and Γ3={x=0,y=0},the following nonlinear operatorFwill be considered:

wherenis the outer unit normal vector.The existence of the Neumann boundary traceKorfollows from the local regularity results in[17]by assumingfm∈L2(Ω).Moreover,to adopt the domainD(F)as well as the observation data spaceY,we confine the nonlinear operatorFmapping from

to

The aim of the current section is to show that the nonlinear operator(2.6)is compact.

The following proposition holds true in view of the estimates(2.3)–(2.4).

Proposition 2.1Let Assumption2.1hold true.Thenanddepend continuously on the exchange rate function α.

ProofDefine two exchange rate functionsα1andα2and by(hm,1,hc,1)and(hm,2,hc,2)denote the solutions of(2.1)–(2.2),respectively.Taking the difference between these two systems and introducingwe obtain

with homogeneous Dirichlet boundary conditions.A basic estimate then holds true in view of(2.3)for the coupled system(2.7)such that

whereCis a constant independent ofα1,α2,fmandfc.

Following the arguments in[13],both integrals on the right-hand side define a bounded linear functionaland have an estimate,by using(2.8),

Since the source termsfmandfcare fixed,we then derive for(2.7),

by implementing the classical elliptic regularity in domains with corners(see[12])and interpolation inequalities(see[1]).The proposition is thus proven.

The next proposition verifies that the local Neumann boundary tracealso depends continuously on the exchange rate functionα.

Proposition 2.2Let Assumption2.1hold true.Thendepend continuously on the exchange rate function α.

ProofDefine two exchange rate functionsα1andα2,the validity of the Neumann boundary data forhm,1andhm,2relies on the local regularity of the CCPF model,i.e.,[17,Lemma 2.3].

Using the same argument as in[10,Section 6.3.1]with respect to the system(2.7),we can prove for an open setandW∩Ωc=∅,an interior estimate such that

The proposition is valid due to the trace theorem.

We then present the following theorem which confirms the compactness of our proposed forward operatorF.

Theorem 2.1Let Assumption2.1hold true.The forward operator F defined in(2.6)is a compact mapping from(L2(Ωc)∩Π)to(L2(Γ1),L2(Γ2),L∞(Γ3)).

ProofThe proofis straight forward by noticing Proposition 2.2 and the compactness of the identify operator from(Hs(Γ1),Hs(Γ2))to(L2(Γ1),L2(Γ2)).

The compactness of the forward operatorFthus indicates the ill-posedness of the inverse problems and calls for the use of regularization methods,for instance iterative regularization schemes,which are presented in the forthcoming section.

3 Parameter Identification Problem for the CCPF Model

3.1 Inverse problems and iterative regularization schemes

The focus of the current section is the identification of the exchange rate functionα(x)defined on Ωcfrom the measurements of the partial Neumann boundary data

i.e.,the boundary seepage velocity(see Figure 3).

These types of parameter identification problems for partial differential equations are typically ill-posed(see[9,14–15]),provided that the forward operatorFis compact.The illposedness means that the solution of(2.5)may not be unique or may not depend continuously on the measurement data,ifit exists.Throughout the current work,we assume that the exact dataz∈Yis attainable,such that there exists an exchange rate functionα†∈(L2(Ωc)∩Π)satisfyingF(α†)=z.In practice,the exact right-hand sidezis not known precisely such that only a noisy observationzγ∈Yis available with

whereγis the noise level.The parameter identification problem then aims at reconstructingαfrom the noisy observation datazγsuch that

holds in an approximate sense.Due to the discontinuity ofF−1,standard algorithms for well-posed problems,e.g.,the Newton iteration,are not suitable for solving(3.2)(see[15]).Regularization methods are then necessary for avoiding an arbitrarily large amplification of the data noise or round-offerrors.

Figure 3 The parameter identification problem(3.2)where the observation data is defined on(Γ1,Γ2,Γ3)and the exchange rate function α to be reconstructed is defined on Ωc.

Iterative regularization schemes are popular in solving nonlinear ill-posed problems.We refer to[15]for a comprehensive discussion of these schemes as well as a well-developed convergence analysis.To incorporate the convex constraint of the subset Π from Assumption 2.1 we define a metric projection such thatPΠ(α)=max(α,ζ).The simplest scheme among all iterative regularization methods is the Landweber iteration which,taking the metric projectionPΠinto account,is read as

An alternative and faster regularization scheme is the Levenberg-Marquardt method.The idea of this method is to find the next iterate as a minimizer of

which gives rise to the Euler equation,by additionally adding the metric projectionPΠ,

whereIis the identity matrix andϵkis a stabilizing parameter.The initial guessα0contains available a priori information about the unknown exact solution.For further reading on the convex constraint and metric projection we recommend[20,Chapter 9]and a recent paper[19].

3.2 Linearization of the forward operator

In both schemes(3.3)–(3.4),one needs to implement the Fr´echet derivativeF′of the forward operatorFas well as its adjoint operator.In the sequel,we provide a formal calculation of these operators for our forward operator(2.6).Consider some test functionv(x)∈(L2(Ωc)∩Π),e.g.,linear B-splines for a representation ofα(x)from the solution space.Then,the linearization of the nonlinear operatorFatαis given by

whereumanducdenote the solutions of the linearization of the CCPF model(2.1)–(2.2)in the directionv,that is,

with homogeneous Dirichlet boundary conditions

The linearized system(3.5)–(3.6)is easily obtained by defining elementsα1(x)=α(x)andα2(x)=α(x)+ηv(x)with a scalarηand by taking the difference of both systems with

Now we verify thatF′(α)is indeed the Fr´echet derivative ofFatα.To start,we define

Standard calculation shows thatsolves the following system,withi=1,2,

and the boundary condition

By the definition of(um,uc)and the appropriate estimate(2.3),we can derive

which yields the following theorem.

Theorem 3.1The forward operator F subject to(2.1)–(2.2)is Fr´echet differentiable.The derivative F′(α)maps v onto the solution of(3.5)–(3.6).

3.3 Adjoint of the linearization forward operator

In an abstract form the underlying CCPF model(2.1)forh=(hm,hc)can be written as

with the operator matrices

Then,the linearized system(3.5)foru=(um,uc)is read as

Moreover,the following proposition holds true.

Proposition 3.1The operator A is self-adjoint.

ProofWe takeh1=(hm,1,hc,1)andh2=(hm,2,hc,2).It follows that

which proves the proposition.

Suppose that we have an elementr=(r1,r2,r3)withr1∈L2(Γ1),r2∈L2(Γ2)andr3∈L∞(Γ3)belonging to the spaceYof our observation(which finally plays the role of the residual in the iteration).In view of Proposition 3.1,we can obtain the adjoint equation of(3.5)–(3.6)satisfying

with the boundary conditions

We note that the adjoint operator ofF′(α)satisfies the variational equality

Notice that the left-hand side in(3.9)yields

Furthermore,by standard calculation,we have

Summarizing,we obtain the variational equality for any test functionvsuch that

The iterative regularization scheme(3.3)or(3.4)together with the variational equality(3.10)allows us to update the iterates by means of the test functionv.

4 Numerical Examples

In this section,we apply both the Landweber iteration(3.3)and the Levenberg-Marquardt method(3.4)for solving the parameter identification problem(3.2)arising in the CCPF model(2.1)–(2.2).

The forward problem solver is realized by the finite element method.Detailed discussions can be found in[5].For convenience,we set the hydraulic conductivityK,the Poiseuille constantDto unity.The exact exchange rate functionα†takes the form of

along the trace Ωc.Such a choice is consistent with the observation that the exchange rate functionαshall be in the order of the hydraulic conductivityK(see[2–3]).In the current work we are interested in a non-constant exchange rate function since the constant one can be approximated by measuring the data near the boundary∂Ωcwhich serves as an initial guess for the iterative regularization schemes.

We adjust the source termsfmandfcin the domain Ωmand Ωcrespectively,so that the forward solution(hm,hc)satisfies

Recall the inverse problem

The noiseξis added to the exact observationzin the following sense by a uniform noise generator with mean zero and=1 such that

where the constantγplays the role of noise level and is set to 10−2.Due to the ill-posedness of the inverse problems,we need an appropriate stopping criteria so that the iterative regularization schemes enforce stability and noise propagation is avoided.Here we simply implement the discrepancy principle(see[15]),i.e.,where the iteration will terminate at stepk∗=k∗(γ,zγ)when the following criterion is satisfied

withτ≥1.In all examples,to start the iterative regularization schemes,we choose the initial guessα0=2 along the whole trace Ωc.The metric projectionPΠis defined with a threshold valueζ=10−16,i.e.,PΠ:=PΠ(·,10−16)in(3.3)–(3.4).The test functionsvfor the exchange rate function are 32 equally distributed linear B-splines.

4.1 Performance of the Landweber iteration

The first example is devoted to the exact observation dataz,where we want to demonstrate the convergence of the proposed Landweber iteration(3.3).In each iteration,the next iterate is updated along test functions following the variational equality(3.10)by substituting the variablerwith the iterative residual

Since there is no noise in the observation data,the stopping criteria is considered when convergence is observed in the sense that

The upper panel in Figure 4 presents the approximated iterative solutionα1584which satisfies the mentioned stopping criteria(4.2).The middle and lower panels in Figure 4 collect the iterative residualas well as the iterative solution errorfork=1,···,1584(x-axis)in the log-scale.As one can observe in these panels,both the iterative residual and the iterative solution error decrease as the iterative stepkincreases.The approximate iterative solutionα1584,which fulfills the stopping criteria,has a solution error≈0.8222.Actually,one can tune the stopping criteria with a smaller value to improve the accuracy but the computational cost increases as well.

Figure 4 Landweber iteration for the reconstruction of the exchange rate function α†with noisefree data.The stopping criteria is(4.2).Upper figure:The approximate iterative solution α1584 versus the exact solution.Middle figure:The iterative residual for k=1,···,1584(x-axis)in the log-scale.Lower figure:The iterative solution error for k=1,···,1584(x-axis)in the log-scale.

We continue with the second numerical test on noisy observation datazγ.The Landweber iteration takes the same form as for the noise-free data but one has to choose the stopping index satisfying the discrepancy principle(4.1).In our implementation,we fix the constantτ=1.01.Figure 5 collects all the numerical results for the noisy data withγ=10−2.Though the approximate iterative solution looks similar to that of the exact data,the iterative residualon the other hand,never breaks the threshold value 10−2(see the middle panel in Figure 5).The discrepancy principle provides an approximate solutionafter 1618 iterative steps with a solution error≈0.7841.

Figure 5 Landweber iteration for the reconstruction of the exchange rate function α†with noisy data γ=10−2.The stopping criteria is the discrepancy principle(4.1).Upper figure:The approximate iterative solutionversus the exact solution.Middle figure:The iterative residualfor k=1,···,1618(x-axis)in the log-scale.Lower figure:The iterative parameter errorfor k=1,···,1618(x-axis)in the log-scale.

The overall computational cost for one Landweber iteration includes calls of the direct and adjoint solvers,i.e.,referring to(3.3)where the forward operatorF(α)and its linearized adjoint formF′(α)∗are necessarily realized.However the number ofiterations can be rather large and the convergence might be slow(see Figure 4).Thus we will consider a faster Newton type method in the coming subsection.

4.2 Performance of the Levenberg-Marquardt method

The Levenberg-Marquardt method(3.4)is computationally realized by solving in each iteration the linear systemin order to obtainbefore applying the projection operator.A numerical discretization of the exchange rate function as,e.g.,a linear combination of linear B-splines then corresponds to a matrix representation of the linear mapping

We test the Levenberg-Marquardt method for the noise-free dataz.The stopping criteria is the same as for the Landweber iteration presented in(4.2)and we collect all the numerical results in Figure 6.As one can see from Figure 6 the residual decreases much faster compared with the Landweber iteration in Figure 4.At the same time,the Levenberg-Marquardt method provides a solution error of≈0.4019.One additional observation is that the approximate iterative solutionα986in the upper panel offigure 6 has a smaller error in the left part of the domain(i.e.,x∈(0,0.5))in comparison to the right part of the domain(i.e.,x∈(0.5,1)).This is because we impose the observation data near the left boundary of the domain Ωmwhich influences the left part of the approximate solution more than the right part.

Similar to the previous subsection,the performance of the Levenberg-Marquardt method for noisy datazγis summarized in Figure 7.This time,the iterative regularization method terminates only after 3 steps where the iterative approximate solution is presented in the upper panel offigure 7.The solution error is≈0.7358 which is better than that of the Landweber iteration but with significantly less computational costs.

Figure 6 Levenberg-Marquardt method for the reconstruction of the exchange rate function α†with noise-free data.Captions of the panels are the same as in Figure 4 with k=986.

Figure 7 Levenberg-Marquardt method for the reconstruction of the exchange rate function α†with noisy data γ=10−2.Captions of the panels are the same as in Figure 5 with k=3.

Finally we will illustrate the ill-posedness of the inverse problem as well as the importance of the discrepancy principle as a stopping criterion.The problem setting is the same as in Figure 7.At the same time,we disable the discrepancy principle and let the iteration proceed until 100 iterations.As one can observe in Figure 8,though the iterative residualmonotonically decreases as the iteration index increases,the solution errordoes not proceed in the same manner.

Figure 8 Levenberg-Marquardt method for the reconstruction of the exchange rate function α†with noisy data γ=10−2but without application of the discrepancy principle.Captions of the panels are the same as in Figure 5 with k=100.

4.3 Further discussion of the exact exchange rate function

Figure 9 Levenberg-Marquardt method for the reconstruction of the exchange rate function α†with projection PΠ.Captions of the panels are the same as in Figure 6 with k=2038.

Recently the authors of[6]revealed that to fit the modified CCPF model with the Stoke-Darcy system by minimizing the difference of the solutions via two models one should choose the nearly optimal choice of exchange rate function sufficiently larger than 25.We also tested an additionalα†(x)=26+10x(1−x)in a similar process.The numerical results are consistent with those in previous subsections and are omitted here.

In both tested cases the metric projectionPΠin(3.3)–(3.4)did not play a role since both exact exchange rate functionsα†are quite large.To validate such a projection operator,we choose a smallα†=0.99∗10−5−0.9∗10−5sin(πx).In the first 6 iterations of the Levenberg-Marquardt method,the operatorPΠprojects the exchange rate functionto the subset Π withζ=10−16.It is observed that the algorithm still performs well and the approximate solution is presented in the upper panel offigure 9 with a solution error≈3.1199∗10−5.For the case of noisy data,situations are similar and we omit the details.

[1]Adams,R.A.,Sobolev Spaces,Pure and Applied Mathematics,Vol.65,Academic Press,New York,London,1975.

[2]Bauer,S.,Liedl,R.and Sauter,M.,Modelling of karst development considering conduit-matrix exchange flow,calibration and reliability in groundwater modelling,Proceeding of the ModelCARE 99 Conference,2000,10–15.

[3]Bauer,S.,Liedl,R.and Sauter,M.,Modeling of karst aquifer genesis:In fluence of exchange flow,Water Resour.Res.,39,2003,1285,12 pages.

[4]Birk,S.,Liedl,R.,Sauter,M.,and Teutsch,G.,Hydraulic boundary conditions as a controlling factor in karst genesis,Water Resour.Res.,39,2003,1004,14 pages.

[5]Cao,Y.Z.,Gunzburger,M.,Hua,F.and Wang,X.M.,Analysis and finite element approximation of a coupled,continuum pipe- flow/darcy model for flow in porous media with embedded conduits,Numer.Methods PDE,27,2011,1242–1252.

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

[7]Chen,N.,Gunzburger,M.and Wang,X.M.,Asymptotic analysis of the differences between the Stokes-Darcy system with different interface conditions and the Stokes-Brinkman system,J.Math.Anal.Appl.,368,2010,658–676.

[8]Discacciati,M.and Quarteroni,A.,Navier-Stokes/Darcy coupling:Modeling,analysis,and numerical approximation,Rev.Mat.Complut.,22,2009,315–425.

[9]Engl,H.W.,Hanke,M.and Neubauer,A.,Regularization ofinverse Problems,Kluwer Academic Publishers,Dordrecht,1996.

[10]Evans,L.C.,Partial Di fferential Equations,Providence,RI:American Math.Soc.,1998.

[11]Faulkner,J.,Hu,B.,Kish,S.and Hua,F.,Laboratory analog and numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains,J.Contam.Hydrol.,110,2009,34–44.

[12]Grisvard,P.,Elliptic problems in nonsmooth domains,Vol.24 of Monographs and Studies in Mathematics,Pitman Advanced Pub.Program,Boston,1985.

[13]Hua,F.,Modeling,Analysis and Simulation of Stokes-Darcy System with Beavers-Joseph Interface Condition,Ph.D.Thesis,Florida State University,Tallahassee,2009.

[14]Isakov,V.,Inverse Problems for Partial Di fferential Equations,2nd ed.,Springer-Verlag,Heidelberg,2006.

[15]Kaltenbacher,B.,Neubauer,A.and Scherzer,O.,Iterative Regularization Methods for Nonlinear Ill-Posed Problems,Walter de Gruyter GmbH&Co.KG,Berlin,2008.

[16]Kiraly,L.,Modelling karst aquifers by the combined discrete channel and continuum approach,Bulletin d’Hydrog´eologie,16,1998,77–98.

[17]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,2012,065003,16 pages.

[18]Shoemaker,W.,Kuniansky,E.,Birk,S.,et al.,Documentation of a conduit flow process(CFP)for MODFLOW,U.S.Geological Survey Techniques and Methods,6.2005,A24.

[19]StÜck,R.,Burger,M.and Hohage,T.,The iteratively regularized Gauss-Newton method with convex constraints and application in 4Pi-Microscopy,Inverse Problems,28,2012,015012,16 pages.

[20]Vogel,C.R.,Computational Methods for Inverse Problems,SIAM,Philadelphia,2002.

[21]Wang,X.M.,On the coupled continuum pipe flow model(CCPF)for flows in karst aquifer,Discret.and Contin.Dyn.Syst.Ser.B,13,2010,489–501.