12
1,
2
,,,
nn
h ijkhijk
n
hijk
nnn
zh h
zijk
ijk
CC
Ct
hUDD CC
kijk kij



 



11
,,,
nn
h ijkijkijkh
CgX

(10a)
(10b)

where


11,
,
nijk ijk
Ux
c U11
1,1,1/ 2,1,1,1/ 2,
1, 1,
1,
nn n
ijk ijkijk ijk
ijkijkxijk
x
HUDDHUD Dc









2
11
2,2, ,1/2,2,2, ,1/2,
2, 2,2,
,1,
nnn n
ijk ijkijk ijk
ijkijkijkijky ijk
y
Uy
cUHUD DHUDDc









3
11
3, 3,,1/23,3,,1/2
3, 3,3,
,1,
nnn n
ijkij kijkij k
ijkijkijkijkz ijk
z
Uz
c UHUDDHUDDc







and

1,0, .
0, 0.
z
Hz z



12/3
,,
,
1
3
12
,
,,,
nn
hijk hijk
s ijk
nnn
hz hh
zijk
HH
St
KCHH
kijk kij




11
,,.
nn
h ijkijkijkh
HhX

Next, for fluid Equation (1) the fractional steps sinite
difference scheme is given by









,
,,
0
 
1/3
,, 1/ 3
,1
23
1
,, ,,
nn
h ijkh ijknn
sijkhx h
xijk
nn nn
hyh hzh
yz
ijk ijk
n
nn h ijk
h ijkh ijkn
ijk ijk
nn
HH
SKCH
t
KCHKC H
C
CC q
t
KCCijkii jk







31
2hh
zijk
1/3 1
,,,
nn
h ijkijkijkh
HhX


(11a)
(11b)
 

 
2/3 1/3
,,
,
2/3
2
12
,
,,,
nn
h ijkhijk
s ijk
nn n
hy hh
ijk
HH
St
KC HH
kjjik




2/3 1
,,,
nn
h ijkijkijkh
HhX


y
ji
(12a)
(12b)
(13a)
(13b)

The initial approximation reads

00
,0, 0
,,
h ijkijkh ijkijkijkh
CcXHHXX .
,
nn
CH n
t
(14)
The algorithm for a time step is as follows: Assuming
the Approximate solution
,,hi
jkhijk at time is
known, it is needed to find out the approximate solutions
at the next time level. First, compute Darcy’s velocity
n
U, from schemes (8a), (8b), and method of speedup is
used to get the solution of transition sheaf
1/3n
C

2/3n
C
,hijk
along x direction. Secondly, from schemes (9a), (9b) we
obtain the solution of transition sheaf ,hijk . Thirdly,
from schemes (10a), (10b) we obtain solution
1n
C
1/3n
H
,hijk .
Next, from (11a), (11b), by using method of speedup, we
get the solution of transition sheaf
,hijk along x
direction; from (12a), (12b) we obtain the solution of
transition sheaf
2/3n
H
,hijk . Finally, from schemes (13a),
(13b) we obtain solution
1
,
n
h ijk
H. So a complete time
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
976
step can be taken. At last, it is because of the positive
condition that this finite difference solution exists, being
the sole one.
2.2. The First Order Weighted Upwind Finite
Difference Fractional Steps Scheme
For salt concentration Equation (2), the first order up-
wind finite difference fractional steps scheme is given by






12
1/3
,,
,
23
,,
,,
*,
,,1
,,
nn
nn
h ijkh ijk
n
hijk x
yh z
yZ
ijk
nn
Ux Uy
h ijkh ijk
nnn
ijkh ijkh ijk
CC
t
DCD
CC
qCCi jk

 





3
1/3
1
,
,
2
,,
n
n
x hijk
nn
h
ijk
n
Uz h ijk
CD
C
C
C
i ijk

1,,
ijk h
X
(8a)
1/3
,
nn
h ijkijk
Cg
 (8b)
 

2/3 ,
n
hijk
D
CC

2/3 1,,
nn
ijk h
X


 
2/3 1/3
,,
2
,
12
,,
,
nn
h ijkh ijk
nn
hijkyh
y
CC
C
t
jikj jik



(9a)
,hijk ijk
Cg (9b)
 

12/3
1,
nn
nnn
hh ijk
CC DCC

11
,,,
nn
h ijkijkijkh
CgX


 
,,
3
,
12
,,
,
hijkh ijk
h ijkz
z
C
t
ki
jk kij


(10a)
(10b)
where
 



1,1, 1,1,
1
1
,
nnn n
Ux
ijkijkijkijkx ijk
cUHUHUc
h
 



1, 1, 1,
12
x
nijk i jki jk
Uc c


and 23
,,
,
nn
Uy Uz
ijk ijk
cc

are defined similarly with a con-
stant
0, 1


3
0,1 ,
.
The algorithm is similar to that of Scheme (8)-(13).
3. Convergence Analysis
For brevity we assume 1,hN
T
ijk

,, ,
X
ih,tnt

,nn
WX tWjh khn
,
ijk ijk . Let
hh
π,
H
HcC
 where H and c are the exact
solutions of this problem (1) - (5), and Hh and Ch are the
difference solutions of the schemes (8) - (13). ,, and
1 denote the inner product and the norms on the dis-
crete spaces l2() and h1() corresponding to L2() and
H1() [19,20,29]. First consider the second order scheme.
Theorem I. Suppose that the exact solutions of prob-
lem (1)-(5) satisfy condition:


Adopt the modified method of upwind procedures (8)-
(13). Let the dissectible satisfy relation:
2
tOh .
Then the following error estimates hold:



1, 1,4,
2
,,
,.
Hc WWL W
4, 22 2
,,
H
tctLWHt c


 t LL
 



2
11 2
22
(; )(; );
*2
(;)
hhth
LJhLJh
L
Jl
th
LJl
HHcC dHH
dcCM th


 (15)
where
 
2
2
;;
0
sup ,sup
N
nn
LJS LJS
SS
ntTNtT n
gg gt


g


*
.
H
M
depends on ,
Constant c
1/3n
C2/3n
h
C

and their derivatives.
Proof. First consider the concentration equation. For
Equations (8)-(13), dispel h and , and we
get the following equivalent form:



123
1
11
,, 1
11
,1
1
11
22
2
1
11
33
3
,,,,
,, ,
12
12
12
nnn
nn
h ijkh ijk
nnn
hijkx h
xijk
ijk
nn
yh
yijk
ijk
nn
zh
Zijk
ijk
nnnn
h ijkh ijkhijkijkh
Ux Uy Uz
CC h
CUDDC
t
hUDD C
hUDD C
CCCqC
















 






*,
,
1
1
21
11
1
1
1
22
2
1
11
11
1
1
1
33
3
1
2
2
12
1||
2
12
12
12
nn
ijkh ijk
nn
xh
x
ijk
nn
yt h
y
ijk
nn
xh
x
ijk
nn
zt h
z
ijk
n
ijk
C
h
tUDDC
hUD DdC
hUD D C
hUDDdC
hUD





 
























1
1
1
1
1
33
3
1
1
3
1
1
1
1
1
122
2
1
1
133
3
12
1||
2
12
1,
2
1,,
n
yh
y
nn
zt h
z
ijk
n
x
ijk
nn
xh y
y
nn n
hzth
z
ijk
DC
hUDDdC
h
tUD
h
DCUDD
h
CUDDdC
ijkN


 






 









1,
(16a)
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
Copyright © 2012 SciRes. IJG
977
11
,
,
ijk ijkh
X


1n
.
(16b)
nn
h ijk
Cg
From Equation (2) (tt) and (16) we have the concentration error equations.


 
1
11 22
11
1
,, 1
11
,1
11
11
11
223 3
23
1
,,
,, ,
12
11
22
nn n
nn
hijk hijk
nnn
hijkx h
xijk
ijk
nnnn
yh zh
yZ
ijk ijk
ijk ijk
nn n
h ijkijkh ijk
Uxu xUyu
h
CUDD
t
hh
UD DUDD
Cc C



 







 

 






11
33
11
,
,,,
11
11
11
111
11
11
11
11
222
22
11
22
11
22
12
nnn
nnn
ijkh ijkijk
yUzuz
nnn
xh
xijk
ijk ijk
nnn
yh
yijk
ijk ijk
cCc
hh
uDUD DC
hh
uDUD DC
h


















 







 







 




11
11
11
333
33
1*, 1*,
,,,
11
11
21111
11 22
12
12
11
22
nnn
zh
zijk
ijk ijk
nnnnnn
ijkh ijkijkijkh ijkh ijk
nnn
xy
xy
ijk
h
uDUD DC
qC cqCC
hh
tuDDcuDDd

 

















 
 

 

 

n
t
c




11
11
1
11 22
12
11
11
111
11 22
12
11
22
11
22
nn
nn
xh yth
xy
ijk ijk
nnn
xy
xy
ijk ijk
hh
UDDCUDDdC
hh
uDDcUD DdC

 











 

 


 




 

 

 


 



ijk
n
th

 


11
11
31111
11 22
12
11
11
133 1
31
1
1
1
12
2
{1( (1(
22
11
22
1||
2
nnn
xy
xy
ijk
nnn
zt
zx
ijk
ijk
nn
xh
hh
tuDDc uDD
hh
uD DdcUD
h
DCUD
 
 
 






11n
c

 
 
 
 
 
 
 
 







1
1
1
23
3
1
31,
1||
2
,1,,1,
nn
yh
y
nn
zt hijk
zijk
h
DC UD
DdC ijkN

 




(17a)
10, ,
n
ijkijk h
X

(17b)
where

12
.
1,
nijk
th

1/3n
h
H2/3n
h
H
Next, consider the fluid equation. For Equations
(11)-(13), dispel and , and we get the
following equivalent form:
Y. R. YUAN ET AL.
978






 


 


 
1
,, 11
,12
1
,
,,
3
0
21
12
1
13
nn
h ijkhijknn nn
sijkhxhhyh
xy
ijk ijk
n
nn h ijk
h ijkh ijknnn
ijkijkh h
zijk
nnn
hxs hyth
xy ijk
nn
hxs hzth
xz
HH
SKCHKCH
t
C
CC qKCC
t
tKCSKCdH
KCSKC dH
 
 
 



 



1
3
nn
hzh
z
ijk
KCH








1
23
311
123 ,1
nn
hys h
yz
ijk
nnnn
hxshys hzth
xyz ijk
KCSKC
tKCSKCSKCdH







, ,1,
nn
zthijk
dH
ijkN
 
11
,.
nn
ijk h
HhX


1n
(18a)
,hijk ijk (18b)
From Equation (1) (tt) and (18) we have the fluid error equations.






 

 




1
1
,
11 1
00
11
33
nn
nn ijkh ijk
ijk ijk
nnnn
hhhijkijk
ijk
nn nn
hh
zz
ijk ijk
cC
Kc KCHq
t
Kc cKCC


 




   







1
111
,1 2 3
ππ
πππ
nn
ijk ijknnnn nn
sijkhxhy hz
xy z
ijkijk ijk
n
ijk
SKCKCKC
t
q
 






 


 


 


 


 


 
2111
12
11
12 1
11
13 2
1
23
nnn
xs yt
xy ijk
nnn n
hxs hythxs
xyx z
ijk
nnn n
hxs hzthys
xzy z
ijk
nn
hys hz
yz
tKcSKcdH
KCSKC dHKcSKc
KCSKC dHKcSKc
KCSKC




 













11
3
11
3
nn
zt ijk
nn
zt ijk
dH
dH


 
The Theory and Application of Upwind Finite Difference Fractional Steps Procedure for Seawater Intrusion
International Journal of Geosciences, 2012, 3, 972-991
http://dx.doi.org/10.4236/ijg.2012.325098 Published Online October 2012 (http://www.SciRP.org/journal/ijg)
The Theory and Application of Upwind Finite Difference
Fractional Steps Procedure for Seawater Intrusion
Yirang Yuan, Hongxing Rui, Dong Liang, Changfeng Li
Institute of Mathematics, Shandong University, Jinan, China
Email: yryuan@sdu.edu.cn
Received August 14, 2012; revised September 10, 2012; accepted October 12, 2012
ABSTRACT
Numerical simulation and theoretical analysis of seawater intrusion is the mathematical basis for modern environmental
science. Its mathematical model is the nonlinear coupled system of partial differential equations with initial-boundary
problems. For a generic case of a three-dimensional bounded region, two kinds of finite difference fractional steps pro-
cedures are put forward. Optimal order estimates in norm are derived for the error in the approximation solution.
The present method has been successfully used in predicting the consequences of seawater intrusion and protection
projects.
2
l
2
l
Keywords: Seawater Intrusion; Three-Dimensional Region; Upwind Fractional Steps; Norm Estimate; Numerical
Simulation
1. Introduction
Seawater intrusion refers to the invasion of salt water
into the groundwater in coastal areas caused by the
changes in natural water environment and social and
economic development. In recent years, it has occurred
in many countries in the world such as USA, Netherlands,
Israel and Japan. After 1970s, the northern coastal area of
our country, especially economic zones around Bohai
such as Shandong, Hebei and Liaoning, is getting more
and more seriously affected by this problem with Shan-
dong province standing out. It leads to the great decrease
in industrial and agricultural production, making people’s
living conditions, especially their drinking water, poorer
and poorer. Therefore, it is urgent that seawater intrusion
be completely tackled.
The mathematical model consists of water head equa-
tion and salt concentration equation. Because of the
compressibility of porous media and that of the fluid, the
change in fluid density with the salt concentration, and
with the consideration of the fact that the salt is in the
moving state in porous media, there may occur convec-
tion-dominated diffusion. While water is moving in the
water-bearing stratum, it carries salt. The movement of
this solute with underground water is called solute con-
vection. Since salt is inhomogeniously distributed in the
whole solution, it always diffuses from places with high-
er concentration to places with lower concentration even
if the solution does not move.
1.1. Water Head Equation
With Darcy’s law, Euler method and Huyakorn’s lin-
earization method, the water head equation is obtained
[1,2].


 
3
0
,,,,0,,
s
T
H
SKHc
tcqxyztJT
t

 
 
e
(1)
where
s
S

is the specific retention,
H
0
pgz

is
water head function, p stands for pressure, 0
repre-
sents the density of reference water (fresh water), g is
gravitational acceleration, z is the height of water con-
taining layer,
is density and depends only on the
concentration of salt c, Hugakorn’s linearization
1
0
s
cc
 
 is adopted. cs is the concentration
corresponding to the maximum density, and
is the
density difference ration
00s


. Kg
,
is viscosity of the fluid,
s
c
is the density
coupling coefficient. e3 = (0,0,1)T,
is the porosity, q
is source or sink term, and the permeability is noted by
1
2
3
00
00.
00
K
KK
K
Cl
1.2. Salt Concentration Equation
dissolved in the fluid causes The movement of
C
opyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 973
convection and diffusion of in porous media. From
Fick’s law and mass conservation law we have the fol-
lowing concentration Equation [1,2].
Cl


T
,,,, ,
Dc c
cxyzt J

0
*
c
t
qC

 

u
Cl*
C

(2)
where c stands for the concentration of , is the
salt concentration near the source well. Darcy’s velocity
and the diffusion matrix are denoted by

3
0
K
Hc
 ue
12 13
22 23
32 33
.
DDD
DD
DDD







, ,,,0
, ,,.

and
11
21
31
DD
1.3. Initial Boundary Value Conditions
To make a complete system, appropriate initial boundary
value conditions are necessary in addition to the above
equations. The initial value condition is

,0

0
0
,, ,,
,,
H
xyH xyz
cxyz
z

,t
cxyz
xyz

1
,,
, .
yzt
(3)
There are three kinds of boundary value conditions.
When concentration and water head are known, the first
kind of boundary value condition can be given as


,,,,, ,,
,,, ,,,
H xyzhxyztcx
g
xyztxyzt J

2,ztJ 

0, Dc
(4a)
The second kind of boundary value condition can be
given to non-flow boundary:
0,, ,xy un n (4b)
where n is the unit vector in outer normal direction. A
kind of Stefan boundary condition is for free surface
boundary.
The boundary condition of water heat equation:

Hz
 
3
3
,,,
10,, ,,.
xyz
H
H
zxyztJ
t

 
uw

',c c wn
Cl
2
l


,,
Dc
xyz

(4c)
The boundary condition of salt concentration equation:

0
3
1
, .t J

 
 
nu
where w is the permeated fluid flow rate in a unit area
and c is the concentration of in permeated fluid.
In the study of seawater intrusion miscible model, for
the miscible fluid Henry suggested an analytic solution
under the simplified boundary condition and with the
steady-state flow in the homogeneous medium [3]. Segol,
Pinder, Grug, Heinrich and others studied the two-
dimensional cut plane problem [4,5], and Huyakorn,
Gupta and Yapa studied the solving process of the three-
dimensional problem [5,6]. However, their calculations
are made in specifically assumed conditions; therefore,
they can not truly reflect seawater intrusion.
For plane incompressible two-phase displacement
which is assumed to be -periodic, Jim Douglas,
Ewing, Russell, and others have published famous papers
on the characteristic finite difference method and finite
element method to solve the convection-dominated dif-
fusion problems with finite difference method, and to
overcome oscillation and faults likely to occur in the
traditional methods [7-12]. For compressible two-phase
displacement problem, Douglas and others have con-
tributed much to the mathematical model of “infinitesi-
mal compressibility”, numerical method and analysis
[13-16]. Douglas and Yuan discarded periodic conditions,
put forward a new modified characteristic finite dif-
ference method and finite element method, and obtained
optimal order estimation in norm [17-20]. Special
treatment is needed for characteristic lines because the
method of characteristics asks for interpolation and they
may go through the boundaries near the solution regions.
It is necessary to find out the intersection point of cha-
racteristic lines and mesh boundaries and calculate their
corresponding functional values. While such calculation
is designed, we must determine whether characteristic
lines really go through the boundaries in order to decide
whether time steps lengths should be changed. In a word,
the practical calculation is quite complicated [19,20].
For the convection-dominated diffusion problems,
Axelsson, Ewing, Lazarov and others proposed upwind
finite difference method [21-23] to overcome oscillation
and to avoid computation complexity of the characteristic
differential method near boundary meshes. Douglas and
Peaceman applied the alternating-direction method to
numerical reservoir simulation [24,25]. By using Fourier
analysis, they succeeded in proving the stability and con-
vergence according to the constant coefficient [26,27].
This paper, starting from the actual case, puts forward the
modified method of upwind with finite difference frac-
tional steps procedure for seawater intrusion. It can over-
come oscillation and diffusion and computative comp-
lexity. At the same time it can convert a three-dimen-
sional problem into three successive one-dimensional
problems, reducing computation complexity and making
computation practical. Moreover, it increases the space
calculation accuracy to the second order. Some tech-
niques, such as calculus of variations, energy method,
operator-splitting, upwind method, commutative law of
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
IJG
974
2
l
2. The Upwind Finite Difference Fractional
Steps Procedure
multiplication of difference operators, decomposition of
high order difference operators, the theory of prior
estimates and techniques are adopted. Optimal order
estimates in norm are derived to determine the error
in the approximate solution. Thus the difficult problem
has been solved [16,28].
For brevity we consider only the first kind of boundary
value problem and the diffusion matrix
,,Dxyz
of
diagonal form. In order to get the solution by using finite
difference method we use mesh region h instead of
region
Generally, this is a positive definite problem:
Copyright © 2012 SciRes.
 

*2 2
*
**
0,
0,,,,0
ii
KKcK Kct
DDxyztD
 
 

*
**
,1,2,3,
,,,
Ki
xyz

(5)
. On space (x, y, z), let step lengths be h3, xi =
ih1, yj = jh2, zk = kh3,
where *
K
, *
K
, , are constants.
*
Our assumptions on the regularity of the solution of
(1)-(5) are denoted collectively by
*


1, 1,
2
,
,.
W W
LL


4,
222
,HcL W
Ht ct


In this paper
M
and
express general positive
constant and general positive small constant respectively,
and have different meanings.


 
 
12
12
12
,,
,,,, .
,,
hijk
ijk iijk

x
yzjikj jik
kijk kij
 


Let h
h
n
ijk
W
stands for the boundaries of , Xijk = (ih1,
jh2,kh3)T, tn=nΔt, W(Xijk, tn) = . Let
 
 
 
,1,,1,
1/2,
,,1,,,1,
,1/2,
,,1,,1
,1/2
,,2,
,,2
,,2,
nnn
hijkh ijkijkh ijk
ijk
nnn
hijkhijki jkhi jk
ij k
nnn
hijkh ijkij kh ijk
ij k
KCKX CKXC
KCKX CKXC
KCKX CKXC






 ,





 
1211 11
1,1,,,,1,
1/2, ,
nnnnnn n
xhxhhh ijkh ijkh ijkhijk
ijk
ijk
KCHh KCHHHH







(6a)
 
11 11
2,,1,,,,1
,1/2 ,
nnnn n
hyhh ijkh ijkh ijkh ijk
yijk
KCKCHHH H

 
 





12
()
nn
h ijk
H h

(6b)
 
1111
3,,1,,,,1
,1/2 ,
nnnn n
zhzhhijkh ijkh ijkh ijk
ij k
KCKCHHH H

 






12nn
hijk
H h
 (6c)




1111
.
n nn nnn
hxhhyhhzh
xyz
ijkijk ijk
KC HKCHKC H
 

  
nn
hhh
ijk
KC H
,
n
h ijk
H and ,
n
h ijk
C be

,n
ijk
(7)
2.1. The Second Order Upwind Finite Difference
Fractional Steps Scheme
Let the finite difference solu-
tions of
n
t
n n
h
,
ijk
cX
H
Xt and , respectively. If
the finite difference solutions h and C
H
are known,
we find the finite difference solutions C, First, compute the approximate Darcy’s velocity
as follows:
1n
h1n
h
H
at
tn+1.

123
,, T
n nnn
U UUU

 

11
,1, ,, ,1,
11
1/2, 1/2,
,
nn
nn nn
hh
h ijkh ijkh ijkh ijk
nn
hh
ijk ijk
KC
HH HH
hh
CC



0
12
nKC
U











 

22
,, 1,,
0
2
,1/2, ,1
2
nn
nn
hh
h ijkh ijk
n
nn
hh
ij kij
KC KC
HH
UCC


,,,1,
22
/2,
,
nn
h ijkh ijk
k
HH
hh


 
 
 
 

 

 

33
,, 1,
0
32
nn
nn
hh
h ijkh ijk
n
nn
hh
KC KC
HH
UCC

,,, 1
33
, 1/2, 1/2
.
nn
h ijkh ijk
ij kij k
HH
hh



 
 
 
 

 
For salt concentration Equation (2), the modified method of upwind with finite difference fractional steps
Y. R. YUAN ET AL. 975
scheme is given by




3
1/3
,
,
n
nn
xh ijk
nn
yh
ijk
nn
zh
ijk
n
h ijk
z
C
C
C
C

2
,,iijk
1,,
ijk h
X
12
1/3
,,
,
1
1
111
1
1
1
222
2
1
1
333
3
,,
,,
*,
,
12
12
12
nn
nn
h ijkh ijk
n
hijk
x
ijk
y
ijk
Z
ijk
nn
h ijkh ijk
Ux UyU
nn
ijkh ijk
CC
Ct
hUD D
hUDD
hUD D
CC
qC C




















,1
,,
n
h ijkijk
(8a)
1/3
,
nn
h ijkijk
Cg
 (8b)
where

0i
CD

,,1,2,3.
i
Di



2/3 1/3
,,
,
1
12/3
222
2
12
1,
2
(, )(, ),
nn
h ijkh ijk
n
hijk
nnn
yh h
yijk
ijk
CC
Ct
hUDD CC
jikj jik



 



2/3 1
,,,
nn
h ijkijkijkh
CgX

(9a)
(9b)




 
12/3
,,
,
1
11
333
3
cb h46 y55 ff5 fs23 fc0 sc0 ls238">

 


31111 1
123
11
123 2
,
1,,1,
n
th ijk
nnnn
xsys zt
xyz ijk
nnnnn
hxshys hzthijk
xyz ijk
dH
tKcSKcSKcdH
KCSKCSKC dH
ijk N

 
 


 






1
,
1
π0, ,
n
ijkijk h
X

(19a)
(19b)
where

12
2, .
nijk
th

We shall introduce the induction hypothesis:


1, 1,
,0,,0,
nn ht


0
sup maxπ
nL
(20)
where
22 2
1, 0,0,
ππ π,.
nn n
h
 

11
πππ
nnn
t ijkijkijk


We consider fluid error Equation (19). Test error Equa-
tion (19) against and summing it up
y parts, we have b
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 979
 
 
 

 
11 11
12
11
312
11
1
1
00
1
π,ππ,ππ,π
2
()π,ππ,ππ,
,π,π
,π
nnn nnn nn
stthxxh yy
nnnn nnnnn
hzzhxxhyyh
nnnn nn
hhhttt
nn
h
nnn
t
SddtK CKC
KCKCKCKC
KcKCHdtd dt
cC
qqdt
 
 
 


 


 





3
ππ,π
n nn
z z








 


 



 



 


11
33
211 1
12
1
12
31111 1
123
11
123
,π
)
,π
nnnn n
hh t
zz
nnn
xs yt
xy
nnnn
hxshytht
xy
nnnn
xsys zt
xyz
nnnn
hxshyshzth
xyz
KccKCC dt
tKcSKcdH
KCSKC dHdt
tKcSKcSKcdH
KCSKCSKC dH






 
 





 




1
2
,π,
nn
t
dt
 π.
n
t
dt
 

(21)
Now we estimate the terms on the right-hand side of
(21).


11
22
,π
π,
nnnn
ht
nn
2
hh
ht
K
cKC
Mtt


 
Hd t
dt


 
(22a)
22
,π
π,
nn
tt
nn
tt
dd t
M
dtdt
 



(22b)


1
1
00
22
2
,π
π,
nn
h
nnn
t
nn
t
cC
qqdt
M
ttd t




(22c)
 




11
33
222
2
,π
π,
nnnn n
hh t
zz
nn n
ht
K
ccKCCdt
M
ttd t
 
 

 
(22d)
For the fifth term on the right-hand side of (21).




 




 


  



31111 1
123
11
123
31
12
11
122
11
,π
π,π
,π
nnnn
xsys zt
xyz
nnnnn
hxshyshzth t
xyz
nnnn
hxshyt t
xy
nnnnn
shytt
xy
x
tKcSKcSKcdH
K
1
hx
n
CSKCSKCdHdt
tKCSKCdd
SKcKCdHd
Kc KC


 
 
 



 





 





KC




(23)
11
2,π.
nnnn
hxsyt t
y
SKc dHd
 



For the first term on the right-hand side of (23), though
the operators

K

are self-co
,
xy
xy
K

njugate,
positive definite and bounded, space region is cubic.
However, their products are generally incommutative.
Noting that
,,,,
xy yxxxy y
yy xxxyyx


we have
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
980

 



 
 
 
3
tK1
12
2
311
12 ,2 ,1
1/ 2,,1/ 2,,1/ 2,1/ 2,
,,1
1
,1 2
1/ 2,,1/ 2,
π π
π()π
π
nnnn
hxshyt t
xy
Nnn nnnn
hhs ijkxytijkhys ijkhxtijk
ijkijkijkijk
ijk
nn
s ijkhxhytij
ijkijk
CSKCdd
tKCKCS dKCSKCd
SKCKCd
 
 





  


 
 


1
21 ,
, 1/2,1/2,
1
2,12
, 1/2,1/2,, 1/2,
13
,1 1/2 ,
π
ππ.
nn nn
khhxs ijkxytijk
ij kijk
nnn
xhys ijkhh
ij kijkijk
nnn
yx sijkhxt ijkyt ijk
ijk
KC KCSd
KCSKCKC
SKCdd h


 



 

(24)
From induction hypothesis (20) we learn that

n
KC
,
1h,
2,
h
KC
n
 
1
11
nn
s h
SKC
are
pression
(24), the positive definite property of 1
12
,,
,
xh
y
KC

bounded. To the first and second terms of ex
s
K
KS
should
be applied, and high-order difference term πn
x
yt ijk
d

should be separated. By using Cauchy inequality to
get
 

,1 ,12
1/2,,1/2,
21 ,
,1/2, 1/2,
ππ
N
t
nn n
ys ijkhhxhytijk
ijk ijk
nn
hhxsi
ij kijk
SK
CSKCKCd
KC KCS





eliminate the terms concerned, we can

2
2
,1/ 2,
πnn
ijkh ijk
KC


31
12 ,
1/ 2,,1/ 2,
,, 1
{nn
hhs ijkxy
ijkijk
ijk
tKCKCSd




1/2, xt
ijkd

11nn
ijks ijk








12
3
132* 3
*
,, 1
22
1
1
ππ
2
π.
N
nn
jkx y tijksx ytijk
ijk
nn
h
dhKStd h
1
22 2
3
32* 3
*
,, 1
1
ππ() ππ
2
N
nn n
xt ytsxytijkh
ijk
M
ddtKSt dhM
 


 

For the third term of (24), we have
t
 

 

(25a)
 




1
2,
12
, 1/2,1/2,, 1/2,
13
,1 1/2 ,
22
1
ππ
ππ.
Nnnn
hysijkhh
ij kijkij k
nnn
yxsijkhxt ijkyt ijk
ijk
nn
hh
C SKCKC
SKCdd h
Mt


 
 
 
(25b)
Sim
3
,, 1
x
ij
k
tK

ilarly, for the other terms, we can obtain

 

 

 




1
12
222
2
11 11
22
,π
,ππ .
nnnnn
ythxshyt ht
xy
nnnnnn
yszth h
yz
dH KCSKCdHd
KcSKcdHMt t

 
 
 

(26)
Now, we consider the sixth term of the right-hand side of (21).
3111
12
nn
xs
xy
tKcSKc

 



 
 



 

1
3
11
123
22222
4 2
3*3 1
*
,π
1πππ .
2
nn
zt
nnnnn
hxshyshzth t
xyz
Nnnnn
sxyztijkhh
cdH
KCSKCSKC dHd
411 11
12
nn
xs ys
xyz
tKcSKcSK
,, 1ijk
K
St dhMtt

 





  

(27)

 

Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 981
For the last term on the right-hand side of (21).


2
2.
ntt
 
From (21)-(28) we can obtain
22
11
2,πππ
nnn n
thh
dtM



 
 (28)
 


2
*11
222222
14
1
ππ,ππ,π
2
ππ π.
nnnnnnn
hhhh
nnnnn n
hh hh
Sdt KCKC
2
st h h
M
t
dhttdt
 


  
(29)
on error equation. Test
err 11n
ijk ijk
Next, consider the concentrati
or Equation (17) against nn
t ijk


 and sum-
ing it up by parts, we have
 
 
11 2
1
1
111
11
1
11
11
11 1
1
,, ,
,,1
2
,
nn n
nnnnnnn
ht txx
nn
nnnn
hth
Ux UxUy
h
CddtDu D
D
Ccd C

 
 


1 11
2231
23
,1 ,1
22
nn nnnn
yy zz
hh
Du
D uD
  



 









 

 







233
11
,,,
11
11
11
111
11
11
11
11
111
22
1
11
3
,,
11 ,
22
11 ,
22
12
nnn
nnnnn
tht
UyUz Uz
nnnn
xxht
nnnn
yyht
n
cdC cdt
hh
uDuDDCd
hh
uDuDDC d
huD

 
 










 



















11
11
11
3
1*,11*,
11
11
311
11 31
13
4
1
1,
2
,
11,
22
12
nnn
zzht
nnn nnnn
t
nnnnn
xx hzztt
n
huDD Cdt
qCCqC Cdt
hh
tuDDCuDDdd
h
tu
 
 






 
 


 



 


 
 

 

 
 


 
 


 
 
 

1
1
1
11
11
11
13233
33
1
1
11,
22
,
nnnn n
xx hyyhzzt
nn
t
D
hh
DCuDDCuDDd
dt








 



 
 




 



 
 





(30)
First, we estimate the second term on the left-hand side of (30).

1 1
1 1
111111
11 11
1
,1 ,1
2
nnnn nnn
x
hh
DuD DuD
 
 
 







(31a)
1 1
122
1
11
11
1
22
,1 .
2
xx x
nnn nn
xxht
h
Du
D Mtdt
 





 


Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
982
Similarly,

1
1
111
22
2
11
11
11 111
2222
22
,1
2
1,1 ,1
22 2
nnnn
yy
nn nnn nn
yyyy
h
DuD
hh
DuDDuDM
 




 










 


 
 


(31b)
22
.
n
ht
tdt
 


1
1
111
33
3
11
11
11 111
3333
33
,1||
2
1,1 ,1
22 2
nnnn
zz
nnnnnnn
zzzz
h
DuD
hh
DuD DuDM
 
 



 










 


 
 


(31c)
22
.
n
ht
tdt
 

Now, we estimate the terms on the right-hand side of (30). In induction hyhesis (20) n
U is bounded, so we have
pot

1
11
22
2
1
11
,,
nn
nnnnn n
hth
Uxu x
CcdtMUu t
 
2
,.
n
t
tdt

 
ilarly,
 (32a)
Sim



11
22 33
11
,, ,,
,,
nn nn
nnn nnn
htht
UyuyUzuz
CcdCcdt
22
2 2
2
22 33.
nn nn nn
ht
M
UuUutt dt


 




(32b)
For the second term on the right-hand side of (30), we have




11
11
11
111
11
11
11 2
11
333
33
11 ,
22
11 ,
22
nnnn
xxht
nnnnnn
zzht
hh
uDUDDCd
hh
uDUDDC dtMuUt
 
 








 







 





(33)
22
.
n
t
t dt


For the third term, we have

22
).
n
tt
t dt


(34)
We consider the fourth term, and we have
1*
,1 1*,2
,(
nnn nnnnn
qCcqCCdtMt


 





11
11
3(1 2
hh
tU


1
11 22
12
11
11
31
1,1/2,2, ,1/2,1,2,
,1, 2,
,, 1
1 ,
2
11
22
nnnnn
xhyt t
xy
Nnnn
ijki jkijkijk
h ijkijkijkx
ijk
DDCUDDdd
hh
tDDC UDUD
 











 


 


 
 
 
 

2
n
d


11
11
1
2, ,1/2,2, ,1/2,2,1,2, ,1/2,
,2, 1,
1
11
1
1,1/2,2,1,
,2, 1,
11
22
11
22
ytijk
nnn n
i jki jkijkijki jk
xh ijkijkijkytijk
nn n
i jkijkijk
yh ijkijkijk
hh
DD C UDUDdD
hh
DCUD UD
 














1
3.
nn
xtijkxyt ijk
ddh
 






Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 983
On positive definite condition:



*
1
1.C

In induction hyhesis (20) n
U is bounde
*1
** *
** *
,0 ,0
i
DD DC

 
pot d, so we have
1
1
0
12
n
hUD b
0,1,2,3,









11
11
31
11 2
12
3
122
23
*2
** 0
{1 1
22
.
2
2
,
nnnnn
xhyt t
xy
nnn
xythtt
hh
tUDDCUDDdd
tDbdMtdd


2
 






 







 
Thus, for the fourth term we can obtain





22
111
11 1
1111
11122
11 2
,
11 1
22 2
nnnnn
yt t
xy
nnnn n
xyt
xy
hh
DDdd
hh h
uDUDDcuDDdc
 
 



 








 
 
 


 

 

 

11
11
31
11
12
11
22
xh
tU
DDCU

  






3
122222
2*21 1
** 0
,
.
2
n
t
nnnnn
xyth h
d
tDbdM t
  


(35a)
Similarly, for the fifth and sixth terms we have








11
11
31
11 33
13
3
1222
2*2
** 0
2222
11
11
22
2
.
nnn
xh zt
xz
nnn
xyt xzt yzt
nnnn
hh
hh
tUDDCUDDd
tDbd dd
Mt
 
 






 

 

 


 


 

,
nn
t
d
(35b)
For the seventh term we have








1
1
1
11 2
12
1
1
1
233
3
4
222222
3*21 1
** 0
22
1,
2
.
2
nn
xx h
nnn n
yy hztt
z
ijk
nnnnn
xyzth h
h
CUD
h
DCUD Ddd
tDbdM t
 
  


 
 








 
(36)
For the last term
1
1
411
n
h
tUDD




2
14
,.Mth  (37)
For error Equation (30), from (31)-(37) we can obtain
2
1
nn n
tt
dtdt
 
 
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
984
11
11
1 1
2
1
2
nn
y
n
x
D



 



211 11
112
12
11
11
11 11
3311
3
1
1
1
22
2
1,1 ,1
22
,1 ,1
22
,1 2
nnnnn
tx xy
nn nnn
zzx
nn
y
hh
dt DuDDu
hh
DuDDuD
h
DuD


 
  


 

  



 
 

 
 








1
1
1
33
3
22222
2
11
,1 2
.
nnn n
yz z
nnnnnnn
hh t
h
DuD
(38)
2
uUtt dt
 
 





 
For fluid error Equation (29), summing over 0nL
M
and noting that 0
π0
, we have
 
 

000
0
22222
2
111
10
ππ,ππ,π
π,π}
th
hhhhh
n
LL
nnnnnn nnnn
hhhhh h
nn
dtKCKC
211
LnL
LL
K
CKC uUtt



 

 


(39)
For the first term on the right-hand side of (39) we have

 
22
1
ππ.
LLL
n
h
n
11
11
,
nnnn n
hhhht
nn
K
CKC EdtM







t

(40)
By
222
4,
n
hh


(41)
w
nnn
uU M
e have

124
.
LLL
tth
22 222
11
11
001
πnL nnn
tth
nnn
dt dtM
 


 

 

(42)
For concentration error Equation (38), summing over 0nL



and noting that 0
0, we can obtain
1
21
11
11
1,1
LnLL
h
dt DuD
 




1
1
0
11
11
111 11 1
2233
23
11
11
00 000 0
112 2
12
22
,1 ,1
22
,1 ,1
22
L
tx x
n
LLL LL L
yyzz
xxyy
hh
Du
D DuD
hh
DuD DuD
  
 


  



 




 


1
11
11
1
11 1
11
1
11
11
1
1
33
3
1,1 1
222
,12
Lnn nn
xx
n
n
nn
z
hh
Du
DuD
hh
h
DuD
 






 
 













1
00 0
33
3
,1 }
2
zz
h
DuD
 



1
22 2
22
,1 1
22
nn n
yy
Du
DuD









122
12
114
3
31
0
1.
2
L
nn nn
zh
n
huDMttht




 


 
 

 


Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
Copyright © 2012 SciRes. IJG
985
Then,

1
21
11 1
11
1
0
11
111 1
223
23
22 2
14
1
0
1,1
22
,1 ,1
22
π.
LnLLL
tx x
n
LL L LL
yyz
Lnn
h
n
h
dtD uD
hh
DuDD
Mttht
 
 
 

 

 


 
 
 
 







11
11
3
L
z
uD


 
(43)
Combining (42) with (43) it yields

24
.tht
22 2222
11 1
11 11
0 0
πππ
L L
nn LLnn
tt
n n
ddt Mt

 
 

 
 

(44)
Applying discrete Gronwall inequality, we have

 


24
.Mth

(45)
It remains to check induction hypothesis (20). First,
for 00
π0
, (20) is correct. If
1olds. From (45) we have
22 22
11
11
0
ππ
LnnLL
tt
n
ddt






0n, since
nL, (20) h
11
1,
LL 1/2
1,
π
M
h

, then induction hypothesis
r 1.nL
For the first order weighted upwind finite difference
fractional steps scheme, we have the following theorem.
Theorem II. Suppose that exact solutions of problem
(1)-(5) satisfy condition:
4,
,

holds fo
 
1, 1,
,Hc WWLW


2
,
4,22 2
,,
H
tctLWHt c

t LL

.
Adopt the first order weighted upwind procedures
(8 mates
hold:
)-(10), (11)-(13). Then the following error esti

 
11 22
(; )(; );
**
()
.
hhth
22
(; )
L
JhLJhLJl
HHcC dHH
dcCM th


  (46)
4. Numerical Simulation Results and
Analysis
heying area of Longkou city as the model area which has
3-dimensional observation grid. This area is on the left
bank of Huangshui River neighboring with Bohai in the
north and Huangshui River in the ea
meters and the width is 700 meters. Its average thickness
is about 17 or 18 meters. In the upper part of the wa-
ter-containing layer there is relatively fine sand, and in
the lower part—coarse sand with gravel which contains
one, two or three layers of mild clay and sludge of dif-
ferent thickness. We decompose this area into three parts
according to the permeability. The section graph and
plane graph are listed in Figures 1 and 2,
The geological parameters are listed in Table 1, where
No., CP, RWS, SY, DD and ICP denote zone number,
coefficients of permeability, rate of water, specific yield,
dispersion degree and infiltration coefficient of precipita-
tion.
Let 20m, 30m, 1m.
xyz
hhh
th
LJ
l
Considering tp
st. Its length is 3000
respectively.
he comlexity of problem, we select Huang-

We compare our
results, real values and the results of others. A represents
ults. The comparison of graphs of water head and
concentration are listed in Figures 3 and 4, respectively.
The section graphs for water and concentration are
in Figures 5 and 6.
olog
) SY DD (m)
the results of Nanjing University [30], and B represents
our res
listed
Table 1. The geical parameters.
CP (m/d) RWS (m1
No. ICP
K
x
xyy
K zz
K
s
S
y
S
L
T
A 17 15
8.0 10–5 0.075 8.3 0.001 0.3
B 103 22
1.2 0–4 0.13 8.3 0.001 0.3
C 7 7
5.0 –5 0.04 0.08 0.0004 0.3
0.11 0.08 0.0004 0.3
1
10
D 63 17
1.0 10–4
Y. R. YUAN ET AL.
986
Figure 1. The graph.
section
Figure 2. The plane graph.
Figure 3. Curves of water level comparison.
Figure 4. Curves of concentration comparison.
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 987
Figure 5. Sptember. ection graph of waterevel computation in Se l
Figur ber.
From the above we can see that the computation re-
sults are exact, and the algorithm given in this paper is
stable and can be used as the algorithm for the simulation
of large-scaled problems.
5. Consequences of Protection Projects and
Applied Modular Form of Project
Adjustment
5.1. Consequences of Projection Projects
The main water conservanc
trusion include projects for water saving, Yellow River
regulation, water retaining and artificial precipitation.
Their ultimate goals are to increase water supplies and to
decrease the extraction of underground water the produc-
tion of human and animal needed, so that the descent of
underground water level will be slowed down and even
underground water be increased. All this is very effective.
Up to now, the protection project results are mainly em-
pirical and qualitative. We have not seen publications
both in China and abroad about the real salt water and
fresh water movement changes after the projects are si-
mulated with computers. There
uantitative and comprehensive predictions of vari- ous
pr
cts of water-saving project on
seawater intrusion. Take the average precipitation amount
in many years (Refer to “Comprehensive Control Plan
Against Seawater Intrusion in Laizhou Bay Area of
Shandong Province”). Simulate water levels and changes
of salt concentration two months after the peak period
(July-August) in the following four conditions: the pre-
sent pumping out, saving water 10%, saving water 20%
and saving water 30%. Water heads and concentrations
of some wells at initial time are listed in Table 2. The
calculation and comparison results are listed in Tables 3
and 4. The predictive sections saving 20% are
From the above we can see the consequences of water
saving projects are remarkable. During raining seasons
the underground water level rises again quickly. In dry
seasons, its descent is slowed down. So the projects slow
down the migration of salt concentration to fresh water
Table 2. The initial values of water head and concentration.
Well No. Water head (m) Salt concentration (mg/L)
e 6. Section graph of Cl concentration in Septem
y works against seawater in-listed in Figures 7 and 8.
are no publications on the
q
ojects. Now we take watersaving project as an example
to discuss the predictive result of the projects.
Scheme: Keep the present precipitation level. Take
into consideration the effe
at water
1-2 –1.01 3667
2-2 –2.20 3000
3-2 –2.77 377
2.87 100
4-2 –3.10 400
5-2 –3.13 98
6-2 –
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
988
Table 3. The effects of water saving projects on water head.
Well number, Water head
Saving water
1-2 2-2 3-2 4-2 5-2 6-2
0 –0.45 –1.75 –2.04 –2.42 –2.32–2.16
10% –0.34 –1.52 –1.80 –2.14 –2.07–1.93
20% –0.23 –1.31 –1.56 –1.87 –1.81–1.71
30% –0.12 –1.10 –1.33 –1.60 –1.55–1.48
Table 4. The effects of water saving projects on salt concen-
tration.
Well number, Water head
Saving water
1-2 2-2 3-2 4-2 5-2 6-2
0 3871 3044 1521 101 98 99
10% 3753 3088 1507 101 98 99
20% 3725 3032 1493 101 98 99
30% 3696 3027 1479 100 98 100
Figure 7. The water head prediction at 20% water saving in
July-August.
Figure 8. The salt concentration prediction at 20% water
saving in July-Augu
areas.
5.2. Predicting the Consequences of
Underground dams and ith the aim
to regulation ur water anrevseaer
intrusion. Underground cut-off walls are built to stop
undercnt aawelod
the dam are l
and are We can ae l
ability of damhepcns
the ing e. ds
op undercurrent and increase water supply, playing the
pitation supply coefficient, playing the role of retaining
andg water. This dams in
sead egiretain and regulate under-
ground water, increase the heof the fresh water
heads in upper vheseseater
intrusid thyi iman fr
invadedgions inweraches. Finally, the dam in
loweres (ne pt seawattrn.
Tidal es uomdwarte
art on the ground and the underground base. Our
both parts are very useful. The
on-the-ground part prevents seawater coming in with
windstorms, while the underground part prevents sea-
water intrusion in common situations because of its small
permeability. The advantages of tidal barrages especially
obvious in of windy period.
There usually are two kinds of anti-percolator, namely:
lower reaches dam in seawater invaded region, and upper
reaches dam in seawater invaded region.
e sea, or in other places where both salt and
esh water move freely. If a large amount of under-
ground water is pumped out in coastal areas, water level
goes down rapidly. When underground water level is
lower than the average tidal level, seawater intrusion
happens. This is because of the continuity between inland
fresh water and seawater. Underground dams reduce
greatly or completely stops the permeability of auto-
chthonous layer. Therefore, cut-off walls can reduce the
possibility of the seawater in lower reaches intruding
trusion thanks to the combined actions of their own and
ater curtain. Underground dams should be located
far from the upstream of tides with the consideration of
ould he built
who
undem.
st.
Underground Dam and Tidal Barrage
Projects
reservoirs are built w
ndergoundd pent wat
urrend seater, since water had is w an
s
safety
ocated
high.
under the ground, so
say th
the stability
epaget the scontro
s is t key oint. Our pratice idicate
followfour ffectsFirstly, undergroun dam
st
role of saving and regulating water. Secondly, they raise
underground water level and increase artificial preci-
supplyin
awater inv
rdly, the upper reache
ed rons
ight
reaches, relie ing t prent wa
on anus plang anportt roleor seawate
re
reach
lo
ear th
re
coast)
s
usiorevener in
barragare usally cpose of to ps: th
p
analysis indicates that
Lower reaches dams should be built on rivers which
empty into th
fr
inland. Moreover, they can retain and regulate the
drainage of underground water. They stop seawater
in
fresh w
tide actions. Otherwise, tidal barrages sh
se upper part is barrage and the lower part is
rground da
As for the upper reaches dams in seawater intruded
areas, since the intrusion has occurred, the underground
walls must be built at the head of these areas with the
aim to retain underground water and increase the fresh
water head from upper reaches and to make seawater
intrusion stable. This is also useful for inland fresh water
areas far from the coast because these dams prevent the
decrease of fresh water amount going into the sea, and
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL. 989
thus prevent the descent of underground water level in
the upper area of intruded areas. The descent of under-
ground water level accelerates seawater going into the
inner part of fresh water areas.
We predicted the effects of the dams on both upper
and lower reaches on seawater intrusion. We chose the
above-mentioned calculation regions. The results are
shown in Table 5. Figure 9 shows the calculated
concentration comparison curve. Where A is the depth of
the lower reaches dam 0 m. water level of the upper
reaches is –1.5 m. B is the depth of lower reaches dam
2m. Water level of the upper reaches is –1 m. C is the
depth of lower reaches dam 4 m. Water level of upper
reaches is –0.5 m. D is the depth of the lower reaches
dam 6 m. Water level of upper reaches is 0 m.
5.3. Applied Modular Form of Project
Adjustment
We should also apply numerical simulation to make
underground water mechanics serve our goal. As for
water supply, we should study how to make the limited
underground water resources exert the most social and
economic benefits, how to limits underground water
level descent within our control and how make water
supply reach the utmost. As for the protection of natural
resources, we should study how to control pollutant
discharge and prevent underground water being polluted,
and how to keep water quality within the permission of
hygienic standards. Here we propose the optimal
Table 5. Cl concentration computation with the effect of
upper reaches and lower reaches dams (after two months).
Observation point
Computation
condition 1 2 3 4 5 6
A 0.103 0.186 2.148 4.019 10.900 15.046
B 0.103 0.184 2.142 3.999 10.678 14.800
C 0.103 0.182 2.136 3.987 10.460 14.531
D 0.103 0.180 2.132 4.006 10.325 14.358
Figure 9. Curves of concentration comparison (two months).
m
r are 4940
m/d and 4227 m/d, respectively. Taking some observed
d well and
applethods, we optimize and study ad-
salt concentration of Case 1 is shown in Figure 10.
n tables, itat the second
pumll acts mre hily n thirst one as
for the level of rved wells. Thus, we can drawhe
follow concn
1) F the fg we
Tab. Adjmf (
m ber
ethod (linear programming) and numerical method. By
their combined efforts the modular form is optimized and
controlled. Namely, we take underground water vari-
ables (water level, discharge, concentration and so on) in
differential equations as the decision variables, use dif-
ference method change them into linear algebraic equa-
tion groups and introduce them into linear programming
model as the constraint conditions.
Now we perform numerical simulation of a real pro-
ject. Assume that there are two pumping wells lying in
some area, whose quantities of pumping wate
well A near pumping wells as a new observe
ying previous m
justed project modes under different cases.
Let it be supposed that the maximum quantity of
pumping water of each well is never more than 5000
m3/d during winter without any rain. Three cased are
considered here to optimize the quantity of each well
with adjusted computation. The first case is that the
groundwater level of observed well doesn’t decrease
(Case 1). The second case is that the increase of the level
is less than 0.1 m (Case 2). And the last case is that the
increase of the level is more than 0.1 m (Case 3). Nu-
merical data under three cases described above are illus-
rated in Table 6. t
With three cases considered above, prediction for
seawater intrusion problems is shown in Table 7, and the
From the data i is easily seen th
ping weffe
obse
oeavthae f
t
inglusios.
orixed pumpin well and observedell, th
le 6usted ode o water saving projectm3/d).
Puping well num
Quantity
1 2 3
Case 1 5000 1840 6840
Case 2 5000 1620 6620
Case 3 5000 2050 7050
Table 7. Changes of salt concentration in soil under differ-
ent conditions (mg/L).
Observation point
Water head
1-1 2-2 3-2 4-2 5-26-2
Case 1 387530431494 101 98 100
Case 2 387030421491 101 98 100
Case 3 392430561526 101 98 99
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
990
Figure 10. The salt concentration of Case 1.
first pumping well can work in the usual form while the
second one should be strictly restricted by applying some
water saving rules nearby.
2) Needs of the second pumping well should be firstly
considered in Yellow River Diversion Project.
If the locations of pumping wells and observation
points are different, the regulation modular forms are
different too. With the establishment of ecological and
environmental control projects, it is possible to get the
timely and accurate observation data about seawater in-
trusion. Therefore, the established control model can do
co
r
c4,
11e
State Education Comm0030422047),
thal Science Foundation of Shandong Province
(Grant No. ZR20AQ012), and dependent Invation
Fon of Song Univ (Grant N10TS
031e authoo thanof. Ewing Rd
Prang Jir their discussion and sugons.
s of Fluids in Porous Media,” A
can Elsevier Publishing Company, Inc., Amsterdam,
1972.
[2] D. Liang and H. Ruhe em
of the Consequence of Predictive Seawater Intrusion and
tion Ps,. Edoceg e
MathematicsfProvince
ao Ocni Pin 19
[3] HenrfeDon Saate
croachment in Coastal Aquifers,” US Geological Surrey
Water Supply Paper, 1964, pp. 1613-1624.
[4] G. Segol, G. F. Pinder and W.G. Gray, “A Galerkin
Element Technique for Calculating the Transient Position
of the Saltwater Front,” Water Resources Research, Vol.
11, No. 2, 1975, pp. 343-347.
doi:10.1029/WR011i002p00343
[5] P. S. Huyakorn, P. S. Anderson and J. W. Mercer, “Salt-
water Intrusion in Aquifers: Development and Testing of
a Three Dimensional Finite Element Model,” Water Re-
sources Research, Vol. 23, No. 2, 1987, pp. 293-312.
doi:10.1029/WR023i002p00293
[6] A. D. Gupta and D. D. Yapa, “Saltwater Encroachment in
an Aquifer: A Case Study,” Water Resources Research
Vo
do
ntrol model calculation for all projects in the entire
region.
6. Acknowledgements
This work was supported by the Major State Basic Re-
search Development Program of China (Grant No.
G19990328), the National Tackling Key Problems Pro-
am (Grant No. 20050200064), the National Natural g
Siences Foundation of China (Grant Nos. 1077112
101244, 10372052), the Doctorate Foundation of th
ission(Grant No.2
e Natur
09In no
undatiohandersityo. 20
). Th
of. Lish
rs want t
ang fo
k Pr. E an
gesti
REFERENCES
[1] J. Bear, “Dynamicmeri-
Y. Yuan,i, “T Mathatics Model
,
Protecroject” In: EJiang, ., Predinof th
2nd Higher
Qingd
Research o
versity
Shandong
gdao,
94, pp. 1-5. ean Uress, Q
H. R. y, “Efcts of ispersi on lt Wr En-
,
Finite
,
l. 18, No. 3, 1982, pp. 546-556.
i:10.1029/WR018i003p00546
[7] J. Douglas Jr. and T. F. Russell, “Numerical Methods for
Convection-Dominated Diffusion Problems Based on
Combining the Method of Characteristics with Finite
Element or Finite Difference Procedures,” SIAM Journal
acement in
of Numerical Analysis, Vol. 19, No. 5, 1982, pp. 781-895.
[8] J. Douglas Jr., “Simulation of Miscible Displ
Porous Media by a Modified Method of Characteristic
Procedure,” Lecture Notes in Mathematics 912, Numeri-
cal Analysis, Dundee, 1981.
[9] J. Douglas Jr., “Finite Difference Methods for Two-Phase
Incompressible Flow in Porous Media,” SIAM Journal of
Numerical Analysis, Vol. 20, No. 4, 1983, pp. 681-696.
doi:10.1137/0720046
[10] R. E. Ewing, T. F. Russell and M. F. Wheeler, “Conver-
gence Analysis of an Approximation of Miscible Dis-
by Mixed Finite Elements
of Characteristics,” Computer
placement in Porous Media
and a Modified Method
Methods in Applied Mechanics and Engineering, Vol. 47,
No. 1-2, 1984, pp. 73-92.
doi:10.1016/0045-7825(84)90048-3
[11] A. Bermudez, M. R. Nogueriras and C. Vazuez, “Nu-
merical Analysis of Convection-Diffusion-Reaction Prob-
lems with Higher Order Characteristics/Finite Elements.
Part I: Time Discretization,” SIAM Journal of Numerical
Analysis, Vol. 44, No. 5, 2006, pp. 1829-1853.
doi:10.1137/040612014
[12] A. Bermudez, M. R. Nogueriras and C. Vazuez, “Nu-
merical Analysis of Convection-Diffusion-Reaction Prob-
lems with Higher Order Characteristics/Finite Elements.
Part II: Fully Discretized Scheme and Quadrature Formu-
las,” SIAM Journal of Numerical Analysis, Vol. 44, No. 5,
2006, pp. 1854-1876. doi:10.1137/040615109
[13] J. Douglas Jr. erical Method for a
Model for Cosplacement in Po-
n, “Time Stepping along Characteristics for the
rvoir Simula-
and J. E. Roberts, “Num
mpressible Miscible Di
rous Media,” Mathematics of Computation, Vol. 4, No.
164, 1983, pp. 441-459.
[14] Y. Yua
Finite Element Approximation of Compressible Miscible
Displacement in Porous Media,” Mathematica Numerica
Sinica, Vol. 14, No. 4, 1992, pp. 385-400.
[15] Y. Yuan, “Finite Difference Methods for a Compressible
Miscible Displacement Problem in Porous Media,” Ma-
thematica Numerica Sinica, Vol. 15, No. 1, 1993, pp.
16-28.
[16] R. E. Ewing, “The Mathematics of Rese
tion,” Society for Industrial and Applied Mathematics,
Philadelphia, 1983. doi:10.1137/1.9781611971071
[17] J. Douglas Jr. and Y. Yuan, “Numerical Simulation of
Copyright © 2012 SciRes. IJG
Y. R. YUAN ET AL.
Copyright © 2012 SciRes. IJG
991
edia Based on Combining
h Mixed Finite Element
of Numerical Simula-
d Finite Element Me-
ensional Moving Boun-
ethod for the Solution of
Parabolic Problems on Composite
Immiscible Flow in Porous M
the Method of Characteristics wit
Procedure,” In: M. F. Wheeler, Ed., Numerical Simula-
tion in Oil Recovery, Spring-Verlag, Minnesota, 1986, pp.
119-131.
[18] Y. Yuan, “Characteristic Finite Difference Methods
Moving Boundary Value Problem
for
tion of Oil Deposit,” Science in China, Se ries A, Vol. 37,
No. 12, 1994, pp. 1412-1453.
[19] Y. Yuan, “The Characteristic Mixe
thod and Analysis for Three-Dim
dary Value Problem,” Science in China, Series A, Vol. 39,
No. 3, 1996, pp. 276-288.
[20] Y. Yuan, “Finite Difference Method and Analysis for
Three-Dimensional Semiconductor Device of Heat Con-
duction,” Science in China, Series A, Vol. 39, No. 11,
1999, pp. 1140-1151.
[21] O. Axelsson and I. A. Gustafasson, “Modified Upwind
Scheme for Convective Transport Equations and the Use
of a Conjugate Gradient M
Non-Symmetric Systems of Equations,” IMA Journal of
Applied Mathematics, Vol. 23, No. 3, 1979, pp. 321-337.
[22] R. E. Ewing, R. D Lazarov and A. T. Vassilevski, “Fini
Difference Scheme for
te
Grids with Refinement in Time and Space,” SIAM Jour-
nal on Numerical Analysis, Vol. 31, No. 6, 1994, pp.
1605-1622. doi:10.1137/0731083
[23] R. D. Lazarov, I. D. Mishev and P. S. Vassilevski, “Finite
Volume Method for Convection-Diffusion Problems,”
SIAM Journal on Numerical Analysis, Vol. 33, No. 1,
1996, pp. 31-55. doi:10.1137/0733003
[24] D. W. Peaceman, “Fundamentals of Numerical Reservoir
Simulation,” Elsevier, Amsterdam, 1980.
[25] G. I. Marchuk, “Splitting and Alternating Direction Me-
thod,” In: P. G. Ciarlet and J. L. Lions, Eds., Handbook of
Numerical Analysis, Elsevior Science Publishers BV,
Paris, 1990, pp. 197-460.
doi:10.1016/S1570-8659(05)80035-3
[26] J. Douglas Jr. and J. E. Gunn, “Two Order Correct Dif-
ference Analogues for the Equation of Multidimensional
Heat Flow,” Mathematics of Computation, Vol. 17, No.
81, 1963, pp. 71-80.
doi:10.1090/S0025-5718-1963-0149676-2
[27] J. Douglas Jr. and J. E. Gunn, “A General Formulation of
Alternation Methods. Part 1. Parabolic and Hyperbolic
Problems,” Numerische Mathematik, Vol. 9, No. 5, 1964,
pp. 428-453.
[28] R. E. Ewing, “Mathematical Modeling and Simulation for
Multiphase Flow in Porous Media. In Numerical Treat-
ment of Multiphase Flow in Porous Media,” Lecture
Notes in Physics, Vol. 552, 2000, pp. 43-57.
[29] A. A. Samarski and B. B. Andreev, “Finite Difference
Methods for Elliptic Equation,” Science Press, Beijing,
1984.
[30] Y. Xue and C. Xie, “Study on Joint-Surface of Saltwater
and Freshwater in Seawater Intrusion Problem,” Nanjing
University Press, Nanjing, 1991, pp. 43-94.