Applied Mathematics and Computation 129 (2002) 469–480 www.elsevier.com/locate/amc
On an inverse problem of pressure recovery arising from soil venting facilities M. Slodicka, H. De Schepper
*
Faculty of Applied Sciences, Department of Mathematical Analysis, Ghent University, Galglaan 2, B-9000 Gent, Belgium
Abstract The technique of soil venting is commonly used for remediation of the unsaturated zone of the soil, which is contaminated by gaseous organic pollutants in the pores. We study a steady-state model with a finite number of extraction wells, where the air flow field in the subsurface is described by a linear elliptic boundary value problem, with a nonlocal Neumann boundary condition, accompanied of a Dirichlet side condition, containing unknown parameters. We develop a constructive method for the parameter identification, i.e., for the determination of the unknown pressure values at the active wells from their measured total discharges. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Inverse problem; Data recovery; Soil venting
1. Introduction The infiltration of organic solvents into the subsurface leads to contamination of the soil and of the ground water system. When volatile hydrocarbons have entered the unsaturated zone, specialised firms often use the technique of soil venting for remediation, which consists in removing the air from the pores of the contaminated soil by pumping. Thus, an air flow field is created towards the pumps or active wells in the subsurface. Extracted air is replaced by clean one from the atmosphere, entering the soil through an unsealed part of the
*
Corresponding author. E-mail address:
[email protected] (H. De Schepper).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 0 5 7 - 1
470
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
Fig. 1. Unsealed surface.
surface (a so-called passive well). Consequently, the volatile contaminants in the gas-phase are removed from the soil. As the volatilisation process is very slow, the cleaning procedure may take over months or even years. On the other hand, the time required to develop a steady-state air flow is rather short, see [1]. Moreover, because the low volatilisation rate of the contaminants causes the under-pressure at the active wells to be small, compared to the absolute value of the atmospheric pressure, it offers no advantage to pump strongly in order to keep the total discharge at a reasonable level. When the surface of a contaminated site is open to the atmosphere, the situation for a single extraction well may look as in Fig. 1. This is an inconvenient situation, because only the vicinity of the well will be effectively cleaned. The value of this ‘‘effective radius’’ depends on the specific situation: depth of the domain, soil structure, heterogenity, soil surface, . . . For a detailed account, we may refer to [1–3]. In order to enlarge the effectively cleaned area, the soil surface can be sealed with concrete, asphalt or iron plates. In this way the entrance of clean air from the atmosphere through the surface can be prohibited. The possibly leaky surface can be modelled through a leakage term in a Robin type BC.
2. Problem setting Let X R3 be a bounded open domain 1, with Lipschitz continuous boundary oX, and containing n cylindrical active wells ðn 2 NÞ, with boundaries Cai ði ¼ 1; . . . ; nÞ, and m cylindrical passive wellsSðm 2 NÞ, with boundaries Sm n Cpi ði ¼ 1; . . . ; mÞ. We introduce the notation CA ¼ i¼1 Cai and CP ¼ i¼1 Cpi . The external boundary of X consists of two complementary parts CD and CN , where Dirichlet and Robin type BCs are prescribed, respectively. The situation is schematically shown in Fig. 2. 1 We are considering a physical 3D-model, but the mathematical results obtained remain valid in a 2D-case as well. Models with reduced dimension have some practical relevance in particular situations (cf. [4]).
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
471
Fig. 2. Top view of the domain with active and passive wells.
We suppose that the gas in the pores of the soil obeys the ideal gas law. Denoting the pressure by p and introducing the variable u ¼ p2 , the steadystate air flow, induced by the active wells, is governed by a linear elliptic PDE for u, accompanied of appropriate BCs. This PDE reads r ðAdif ru acon uÞ þ asou u ¼ f
in X:
ð1Þ
Here Adif ðxÞ and acon ðxÞ correspond to the diffusive and convective part of the air flow at x 2 X, respectively. Further, f ðxÞ asou ðxÞu represents the source term, supposed to be a linear function of u. A delicate point in the flow modelling is the choice of appropriate BCs at the active wells. Usually, the discharge (total flux through a well) can be measured, while the distribution of the flux itself along the boundary of an extraction well is unknown. Moreover, from the physical point of view we may assume that a constant, yet unknown under-pressure is induced at the boundary of each active well, by a given discharge. Thus, we are led to a nonlocal Neumann BC, accompanied of a Dirichlet side condition at each boundary Cai ði ¼ 1; . . . ; nÞ, i.e., Z ðAdif ru acon uÞ m ¼ si ðgivenÞ; ð2Þ Cai
u ¼ ci ðunknown constantÞ
on Cai
ð3Þ
for a given set of discharges si 2 R ði ¼ 1; . . . ; nÞ. Here m stands for the outward unit normal vector. Finally, BCs must be supplied at CD , CP and CN as well. They may be taken to be u ¼ gDir
on CD [ CP ;
ðAdif ru acon uÞ m gRob u ¼ gNeu
ð4Þ on CN :
ð5Þ
472
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
The problem (1)–(5) can be interpreted as a problem of parameter identification (inverse problem): the unknown parameters ci (squared pressures at Cai ), entering (3), are sought in terms of the set of given (measured) total fluxes si through Cai ði ¼ 1; . . . ; nÞ. The mathematical model (1)–(5) has been introduced in [4], where also the existence of an exact solution has been proved, under suitable regularity conditions on the data. However, in that paper, no method has been given for the construction of the solution. In fact, to avoid the difficulties connected with the nonlocal BCs (2), the authors have used Dirac functions to model the active wells as point sinks, i.e., the original problem (1)–(5) was replaced by an approximate one. We deal with a new proof of the existence of a solution to the problem (1)– (5), which shows some important advantages. First, the proof also includes the uniqueness of the solution. Secondly, the proof is constructive: it provides a construction method for the solution, starting from the solution of wellposed auxiliary elliptic problems with standard (local) BCs. Finally, the method presented can be extended to the case where the active wells operate in groups, a strategy which is often used when designing soil venting facilities, viz when in each of the disjoint groups the wells work at the same specific under-pressure. The effectiveness of the method is illustrated by a numerical example.
3. Constructive solution strategy Consider the function space on X, V ¼ u 2 W 1;2 ðXÞ j u ¼ 0 on CD [ CP ; where W 1;2 ðXÞ stands for the standard, first-order Sobolev space. To simplify the notations, we denote by Gai ðuÞ the total flux through the boundary Cai of the ith active well, i.e., Z ðAdif ru acon uÞ m; i ¼ 1; . . . ; n: ð6Þ Gai ðuÞ ¼ Cai
Next, we put Au ¼ r ðAdif ru acon uÞ þ asou u:
ð7Þ
Moreover, we introduce the bilinear form að:; :Þ on V V by Z Z Z aðu; wÞ ¼ ðAdif ru þ acon uÞ rw þ asou uw þ gRob uw 8u; w 2 V : X
X
CN
ð8Þ
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
473
We assume that the bilinear form að:; :Þ fulfills the standard ellipticity property Z 2 2 9C0 > 0 : aðu; uÞ P C0 juj1;2;X C0 jruj 8u 2 V : ð9Þ X
Remark 1. We formulate a set of conditions on the data functions, appearing in (8), which ensure the V-ellipticity R 2 9Ad > 0 : X Adif ru ru P Ad juj1;2;X 8u 2 V ; gRob P 0 on CN ; 2 asou P As > kacon k0;1;X =ð2eAd Þ in X for some e 2 ð0; 2Þ: Here, k:k0;1;X stands for the usual L1 ðXÞ-norm. By the Cauchy–Schwarz inequality in L2 ðXÞ and a well-known algebraic inequality, we have Z acon u ru 6 kacon k 0;1;X kuk0;2;X juj1;2;X X
2
6
kacon k0;1;X Ad e 2 2 juj1;2;X ; kuk0;2;X þ 2 2eAd
where k:k0;2;X denotes the L2 ðXÞ-norm. On account of this estimate and the assumptions above, we arrive at ! 2 kacon k0;1;X Ad e 2 2 9C > 0 : aðu; uÞ P Ad juj1;2;X þ As kuk0;2;X 2 2eAd 2
P Ckuk1;2;X
8u 2 V :
First of all, we prove the uniqueness of the weak solution of the BVP considered. Theorem 1 (uniqueness). The BVP (1)–(5) has at most one solution. Proof. Suppose that u1 and u2 are solutions of (1)–(5), and denote u ¼ u1 u2 . Then we have Au ¼ 0
u ¼ c~i ðconstantÞ Gai ðuÞ
ð10Þ
in X;
¼ 0;
on
Cai ;
i ¼ 1; . . . ; n;
i ¼ 1; . . . ; n;
ð11Þ ð12Þ
while moreover, u fulfills homogeneous Dirichlet and Robin BCs on CD [ CP and CN , respectively. Multiplying both sides of (10) by u, integrating over X, applying Green’s theorem and invoking (8) leads to
474
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
aðu; uÞ
n Z X
ðAdif ru þ acon uÞ mu ¼ 0 Cai
i¼1
or, on account of (11) and (12), aðu; uÞ ¼ 0. As u 2 V , we infer from the ellipticity condition (9) that u ¼ 0 a.e. in X, i.e., u1 ¼ u2 . To prepare the constructive proof of the existence theorem, we consider the following auxiliary problems, related to (1)–(5): r ðAdif rv acon vÞ þ asou v ¼ f
in X;
v ¼ gDir on CD [ CP ; ðAdif rv acon vÞ m gRob v ¼ gNeu v¼0
on CN ;
ð13Þ
on CA
and for i ¼ 1; . . . ; n, r ðAdif rzi acon zi Þ þ asou zi ¼ 0 zi ¼ 0
in X;
on CD [ CP ;
ðAdif rzi acon zi Þ m gRob zi ¼ 0 on CN ; zi ¼ 1 on Cai ; zi ¼ 0
ð14Þ
on CA n Cai :
From the theory of linear elliptic equations (cf., e.g. [5]), we know that each of these problems has a unique weak solution under classical regularity conditions for the data. Now, for any a ¼ ða1 ; . . . ; an Þ 2 Rn we define ua ¼ v þ
n X
ð15Þ
ai zi :
i¼1
By linear superposition, the function ua is seen to obey r ðAdif rua acon ua Þ þ asou ua ¼ f ua ¼ gDir on CD [ CP ; ðAdif rua acon ua Þ m gRob ua ¼ gNeu ua ¼ ai
in X; on CN ;
ð16Þ
on Cai ; i ¼ 1; . . . ; n:
The flux of ua through the boundary of the jth active well takes the form Gaj ðua Þ ¼ Gaj ðvÞ þ
n X
ai Gaj ðzi Þ:
i¼1
In order to solve the original problem (1)–(5), we are looking for an n-tuple a, for which
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
Gaj ðua Þ ¼ sj ;
475
j ¼ 1; . . . ; n:
This leads to a linear algebraic system 0 a 10 1 0 a 1 0 1 G1 ðz1 Þ . . . Ga1 ðzn Þ G1 ðvÞ a1 s1 B .. .. CB .. C þ B .. C ¼ B .. C: .. @ . A @ A @ A @ . . A . . . a a a an sn Gn ðz1 Þ . . . Gn ðzn Þ Gn ðvÞ
ð17Þ
Now, we are in a position to prove the existence of a solution of the BVP (1)–(5). Theorem 2. There exists a solution of the BVP (1)–(5). Proof. According to the considerations above, it is sufficient to show the regularity of the matrix G ¼ ðGai ðzj ÞÞi;j entering (17). To this end, suppose that G is singular. Then, one of the columns of G is a linear combination of the other columns. Without loss of generality we may assume that, for some constants k1 ; . . . ; kn1 , ! n1 X a ki zi ¼ 0 for j ¼ 1; . . . ; n: Gj zn i¼1
From this relation and from (14), the function w ¼ zn solve the problem
Pn1 i¼1
ki zi is seen to
r ðAdif rw acon wÞ þ asou w ¼ 0 in X; w ¼ 0 on CD [ CP ; ðAdif rw acon wÞ m gRob w ¼ 0 w ¼ const on Gai ðwÞ ¼ 0 on
Cai ; Cai ;
on CN ;
i ¼ 1; . . . ; n:
Clearly, w ¼ 0 is a solution of this problem. By Theorem 1 this zero solution is the only solution. We conclude that zn ¼
n1 X
ki zi :
ð18Þ
i¼1
Recalling the fact that, by (14), z1 ; . . . ; zn1 vanish on Can , we get zn ¼ 0 on Can . This is contradictory with the BC zn ¼ 1 on Can , see also (14). Hence, the matrix G is regular. Consequently, the algebraic system (17) has a unique solution a. The unique solution of the BVP (1)–(5) is then defined by the relation (15). We emphasize the constructive character of the existence proof. The desired solution of the problem (1)–(5) has been expressed in terms of solutions of
476
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
some auxiliary problems with standard BCs. Moreover, the values a1 ; . . . ; an are nothing else than the, a priori unknown, values of the squared pressure at the respective boundaries Ca1 ; . . . ; Can of the active wells. In practice, the auxiliary problems (13) and (14) are solved by suitable approximation methods, such as finite element methods (FEMs).
4. Extension to wells operating in groups As mentioned above, soil venting is a slow process, which can take over months or years. Due to the small volatilisation rates of the organic contaminants, the under-pressures at all active wells should be low. On the other hand, the engines for extraction wells are nowadays very effective and powerful. Hence, it is useful to design the soil venting facilities in such a way that the active wells are splitted into disjoint groups, each group being connected with a single engine, i.e., the under-pressures at all extraction wells inside one group are equal. We suppose to have k groups ðk < nÞ. Let Wj f1; . . . ; ng denote the set of indices of the active wells belonging to the jth group. The discharge conditions (2) and (3) have to be modified into the form u ¼ aj ¼ const ðunknownÞ X
on
[ i2Wj
Gai ðuÞ
¼ sj 2 R
Cai ; ð19Þ
i2Wj
for j 2 f1; . . . ; kg. Recall that all simultaneously working wells inside one group operate at the same under-pressure, i.e., this group has only one degree of freedom. Thus, for j ¼ 1; . . . ; k, we consider the function X yj ¼ zi ; ð20Þ i2Wj
where zi is the unique solution of problem (14). For any a ¼ ða1 ; . . . ; ak Þ 2 Rk , we introduce the function ua ¼ v þ
k X
aj y j :
ð21Þ
j¼1
Here, again, v denotes the unique solution of the BVP (13). By the superposition principle, the function ua is now seen to obey
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
r ðAdif rua acon ua Þ þ asou ua ¼ f
477
in X;
ua ¼ gDir on CD [ CP ; ðAdif rua acon ua Þ m gRob ua ¼ gNeu [ Cai ; j ¼ 1; . . . ; k: ua ¼ aj on
on CN ;
i2Wj
Analogously as before, the flux of ua through GaWj ðua Þ ¼ GaWj ðvÞ þ
k X
S
i2Wj
Cai reads
ai GaWj ðyi Þ;
ð22Þ
i¼1
where GaWj ðwÞ ¼
X
Gai ðwÞ:
i2Wj
To solve the original, modified problem, we are looking for an a 2 Rk , satisfying GaWj ðua Þ ¼ sj ;
j ¼ 1; . . . ; k:
Hence, we have to solve the linear algebraic system 0 a 10 1 0 a 1 0 1 GW1 ðvÞ GW1 ðy1 Þ . . . GaW1 ðyk Þ a1 s1 B CB .. C B .. C B .. C .. .. .. @ A@ . A þ @ . A ¼ @ . A: . . . a a ak sk GWk ðy1 Þ . . . GWk ðyk Þ GaWk ðvÞ
ð23Þ
In a similar way as for separately operating wells we may show that the matrix ðGWi ðyj ÞÞi;j is regular. The solution a of the system (23) determines the pressure at the boundaries of the extraction wells in each of the k groups, al being the pressure corresponding to the lth group ðl ¼ 1; . . . ; kÞ. Then, the solution ua of the modified problem is the function given by (21).
5. An illustrative numerical example In this section, we present a numerical example, representing a model situation in soil venting. We consider a cylindrical domain in R3 with a vertical axis. We assume that this domain is insulated at all sides (at the bottom by an impervious layer, at the top by an iron plate or concrete, and at the lateral surface by an artificially constructed barrier). Fresh air can enter the domain through one passive well which is open to the atmosphere. When we suppose the geological soil layers to be horizontal and the extraction wells to pump air along the whole of their vertical length, a horizontal air flow field will be created. Vertical integration then yields a reduced
478
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
Fig. 3. Situation.
2D-model (cf. [4]). An example of such a reduced 2D-situation is represented in Fig. 3, where three active wells, with respective boundaries C1 , C2 and C3 , and one passive well with boundary CP are considered. The geometrical data are as follows. The domain is centered at the origin and has radius 2. All wells (active and passive) have radius 0.2. The passive well is also centered at the origin, while the midpoints of the active wells are respectively given by ð1:5; 0Þ, ð1:2; 0:5Þ and ð0; 1:0Þ. We consider the following model problem: Du ¼ 0 u¼0
in X;
on CP ;
ru m ¼ 0 on CN ; u ¼ ci ðunknownÞ on Ci ; i ¼ 1; 2; 3; Z ru m ¼ si i ¼ 1; 2; 3 Ci
for s1 ¼ 0:2, s2 ¼ 0:3 and s3 ¼ 0:25. This problem is solved by the method outlined in Section 3. As the auxiliary problem (13) now is completely homogeneous, we have v ¼ 0. To solve the three auxiliary problems of the type (14), we use a FEM, with piecewise linear polynomials on a triangular mesh. The unstructured triangulation of the domain consists of 21 888 triangles. The isolines for zi , i ¼ 1; 2; 3, are shown in Fig. 4. The corresponding algebraic system of the type (17) is found to be 1 0 10 1 0 0:2 1:77707 0:0220422 0:274207 a1 @ 0:0223199 2:14177 0:166005 A@ a2 A ¼ @ 0:3 A: 0:25 0:287506 0:16705 2:63833 a3
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
479
Fig. 4. Isolines for z1 , z2 and z3 .
Its solution reads a1 ¼ 0:1327387057;
a2 ¼ 0:1506593337;
a3 ¼ 0:1187610405:
The equipotential curves for the corresponding solution ua (15), are shown in Fig. 5. To check the reliability of the method, we consider the following elliptic problem, with standard Dirichlet BCs on Ci , given by the values obtained for ai (i ¼ 1; 2; 3): Dua ¼ 0 in X; ua ¼ 0 on CP ; rua m ¼ 0 on CN ; ua ¼ 0:1327387057 ua ¼ 0:1506593337
on C1 ; on C2 ;
ua ¼ 0:1187610405
on C3 :
480
M. Slodicka, H. De Schepper / Appl. Math. Comput. 129 (2002) 469–480
Fig. 5. Isolines for ua .
When solving this well-posed problem by the FEM mentioned above, the (approximate) values of the corresponding fluxes through the boundaries of the active wells are found to be G1 ðua Þ ¼ 0:2;
G2 ðua Þ ¼ 0:299999;
G3 ðua Þ ¼ 0:25:
These calculated fluxes coincide almost exactly with the prescribed values s1 , s2 and s3 , respectively.
Acknowledgements One of the authors (MS) was financially supported by the VEO-project no. 011 VO 697. The authors thank R. Van Keer, coordinator of this project, for stimulating discussions and for his critical reading of the text.
References [1] D.H. Hiller, Performance characteristics of vapor extraction systems operated in Europe, in: Proceedings of the Symposium on Soil Venting, Houston, Texas, 1991; Environmental Protection Agency EPA/600/R-92/174. [2] D.J. Wilson, Modeling of Insitu Techniques for Treatment of Contaminated Soils: Soil Vapor Extraction, Sparging, and Bioventing, Technomic Publishing Company, Lancester-Basel, 1995. [3] S. Schumacher, M. Slodicka, Effective radius and underpressure of an air extraction well in a heterogeneous porous medium, Transport Porous Media 29 (1997) 323–340. [4] H. Gerke, U. Hornung, Y. Kelanemer, M. Slodicka, S. Schumacher, Optimal Control of Soil Venting: Mathematical Modeling and Applications, Birkh€auser, Basel, 1999. [5] D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, Berlin, 1983.