- Email: [email protected]
Nonlinear simulation of environmental resources systems
114
Mathematicsand Computersin SimulationXXIV
(1982) 114-121
North-Holland PublishingCompany
NONLINEAR
SIMULATION
OF ENVIRONMENTAL
RESOURCES SYSTEMS
V.V. NGUYEN and E.F. WOOD Water Resources
Program,
Department
of Civil Engineering,
Princeton
University,
Princeton,
NJ 08544,
U.S.A.
Nonlinearity, inherent in physical systems, can be simulated using elementary catastrophe theory. In this paper, we present a general framework for the application of this theory to environmental systems and demonstrate its practicality in a series of quantitative models consisting of the eutrophication of small lakes, rainfall infiltration process, geothermal simulation and the multispecies/multiphase flow in oil reservoir. INTRODUCTION The application of elementary catastrophe theory originated by Thorn [l] has evolved to a practical stage where the simulation of environmental systems can benefit from its unorthodox features. The state-of-the-art has been reviewed recently by Stewart [Z] and is well documented in [3].
Generally speaking, natural systems are often nonlinear and consequently exhibit qualitative changes in the structure of the state variables as related external or control parameters reach critical thresholds. Linear models, when employed to simulate these systems, fail to depict In order to propthis characteristic behavior. erly understand such systems, the geometric essence of the state must be brought forth according to conservation laws, constitutive axioms and intrinsic relationships between specific combination of external parameters and the state. When the behavior of the system is governed by the principle of structural stability, i.e., the description of the critical threshold remains invariant under perturbation of the governing energy field, the state-parameter relationships can be classified in terms of elementary catastrophes. For a discussion of these basic forms, we refer the reader to [l] "1. [3]. In this paper we explain how to employ elementary catastrophes to describe or numerically simulate environmental systems. A brief summary of a general framework of application is presented and is followed by a series of models. Two well-known forms are used in this presentation: the cusp and the butterfly catastrophes. Technical details of these models are given elsewhere, see 141, 151, [6], 171, and I&31. 1.
SUMMARY OF THE THEORETICAL FRAMEWORK _-
A mathematical description of a physical process subject to conservation laws and constitutive axioms usually has the form
0378-4754/82/0000-0000/$02.75
N(s;
a,b,c,...) = 0
(i.1)
where s is the state vector and ia,b,c,...} is the set of external or control parameters. The nonlinear operator N(e) consists of the time a (*): the local spatial changes, rate change, 5 such as a/a/ ~a dX
CL = (a 1'"" n lcll = c ui,
’ where x is a position vector,
o1 ), axa n
= ax
al a2 ax . ..ax 12
dt
,
and/or a combination of the latter.
When written as a discrete simulator, the form
*
cln
(1.1) has
c;. 2)
= [A]; + [B]
together with the initial conditions, s is the approximation of the state vector s at nodal ;oi;ts, [A] is the coefficient matrix containing a, ,c,... i and [I?]is the column vector of boundary and previous time level information. When the syste? is nonlinear, [A] will also be a function of s. In order to solve for sAin this situation, the relationship between s and [A] must be given a priori. For example, a; nodal points (31, the relationship between s i and [a,b,c,...) may be written as
P(sj; a,b,c,....) =
0
(1.3)
where P is a polynomial in the given variables. Notice that the time scale for relation (1.3) is assumed to be smaller than that of the evolution equation (1.2). Consequently, we can obtain dsj/dt from (;.3) (at least in principle and Since usually by somq iteration scheme). n s z (l.), (1.3) should be used in parallel to (1.2) in determining the approximate solution of (1.1). When we consider the principle of structural stability to establish (1.3), P(e) is usually best described by one of the seven ele-
0 1982 IMACS/North-Holland
V. V. Nguyen, E.F.
Wood /Nonlinear
simulation of environmental
mentary catastrophes of Thorn [I]. Each parameter a, b, c, . . .. may be regarded as a lumped parameter of other physical parameters of the dynamics under study. Let the normal form of P be
F(Sj; a,b,c,d,e) = 0
(1.4)
where the total number of lump parameters may not exceed five. From (1.4), the set of bifurcation values of {a,b,c,d,e> can be established by defining it to be the values that satisfy (1.5)
Essentially (1.5) denotes the coy?ditions by which the qualitative change in s at the smaller j time scale takes place with respect to change in As the external parameters cross m,d,e1. these thresholds in the control space, the method we use to solve (1.2) at the larger time scale must be adjusted accordingly. In other words, the loci of critical values of {a,b,c,d,e),i.e. the bifurcation set,partition the control space into different characteristic zones within which the corresponding behavior of the state s popsesses different structures. 2.
APPLICATIONS
We illustrate the above framework of utilizing the elementary catastrophe theory in the following analysis of environmental systems. 2.1
E; = F(x;
a,b)
(2.4)
By shortwhere E is a small positive constant. However, term, we mean two to three days periods. at a larger time scale, monthly or seasonally, the trend of growth in P is not steady, viz. dP/dt # 0
(2.5)
Therefore (2.4) may be used to specify the critical combination of chemical, mechanical and thermodynamic factors by which P experiences short-term blooms and crashes. The model of algae dynamics is presented with field data in Figures 1, 2, and 3. Chlorophyll a and ammonium nitrate are used as indicators of P and N ree' spectively.
300-
LAKE
555
SUWER 250
1972
-
3200. z
i
-SAMPLEDULUES
I u
no-
-sLorrcHANsE -II*py)owJGE
100 -
0
Let P denote the b&omass of phytoplankton representing the state s, by using the Monad and the Michaelis-Menten relations, we can show that
0
5
IO YAXIYUM
Figure 1. f2’
115
50.
Summer Algae Dynamics in Nonstratified Lakes 141
P@; fl’
resources systems
f3) = - (P3 + flP2 + f2P + f3) = 0
I5 Am
20
TEYPmmJRE.
25
30
-c
Phytoplankton Biomass as a Multivalued-Function of Air Temperature.
(2.1) where
fl
= fl
f2 = f2
(I,
Ne,
R)
(2.2.3)
(I,
Ne,
R)
(2.2b)
f3 = f3 (Ta, I) and
(2.2c)
I= light influx density, = chemical nutrient level in the lake e environment, R= endogeneous respiration factor, and Ta = air temperature. N
Using the transformation
x = P - f1/3, we can
reduce (2.1) to the cusp catastrophe form, i.e. F(x;
a,b) = - (x3 + a(N e,R)x + b(Ta,I))
3
(2.3)
According to the principle of structural stability [3],we may postulate the short-term dynamics of x by
Figure 2.
Control Space Trajectory for Lake 885.
116
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems where a, b, c are functions of oils.
When the
system conditions are known, i.e. 8,'s are deter1 mined a priori, one can find the critical rainfall rate r as a root of (2.9) at which the c upper layer of the porous medium assumes the state of ixmnediate ponding. This result implies different flow regimes associated with the wetting curve, see Figure 4. Denoting I to be the bulk infiltration volume, i.e., I = sdf where dz is the constant thickness of the moisture front, we obtain a morphological the variation of I in Figure 5.
descriotion of
-140 PLAINFIELD
Figure 3.
2.2
Morphological Representation of the state of P on a Cusp Catastrophe.
Physics of Rainfall Infiltration
[5]
The imbibition of water at the soil surface due to rainfall sometimes exhibits bifurcation behavior known as the Haines jumps. Using the physics of flow in the unsaturated zone, the combined equation of movement of air and water can be obtained according to the generalized Darcy's law. The integral form of this equation over a specific range of moisture content becomes (aHc + rt) 2
= P(s; ~1, ~2, p3, ~4)
REGION
where the state s is the normalized moisture content, H, is the effective capillary drive, r is the constant rainfall rate, t is the "shorttern?' time of precipitation, a is a gravitational characteristic constant, and p.'s are the para1 meters representing soil type, wetness condition, rainfall intensity and rainwater properties. Under commonly encountered conditions, the normal form of (2.6) may be approximated by 3
e3 + (0, +F)x
e5 + (e, +T)
(2.7)
are transformed representations t and (aHc + rt), respectively, and oils are
where x and E of
s
the transformed parameters of pi's_ the bifurcation 4(8
2
+2)3 r
From (1.5),
set is determined by
+ 27(8
4
+$2 r
= 0
(2.8)
which leads r
3
+ ar
2
+br+C=O
OF IMMEDIATE
?ONDING
(2.6) X.
dx Etx=X
SAND
(2.9)
Figure 4.
NDFWALIZED
MOISTURE CONTENT
Unsaturated Flow Regimes for Plainfield Sand.
V. V. Nguyen, E.F. Wood jNonlinear
117
simulation of environmental resources systems and
h,(e) = CAPILLARY DRIVE
’
a0
4%
I[
5
c2
c3
c4
=
“ah
ah (l-$)0 2 m ap
(l-$)P,
-4
2 -
4P
(2.12)
1 1 --c
5
-c
0 -c
6
=
7
k
1
K -$ +(; p,p,SvSJ(hv-h&e "V
k - +' ww
aT Kah (2.13)
'8
(2.14)
= $PV
subscripts m : solid phase, v : steam phase and w f liquid water phase,
h,(,&=SOIL
Figure 5.
h = t
PROPERTIES
Morphological Representation of Bulk Infiltration Volume over Capillary Hysteresis.
sv+s
(2.15)
(SvP,hv + SwP,hw)
w
(2.16)
=1
(2.17)
v = ; (S"U"V" + S"PWVJ 2.3
Geothermal Reservoir Simulation with Nonequilibrium Thermodynamics [7], [8]
The state s of a geothermal flow consists of In one dimension of pressure and enthalpy. space, the flow equations are represented by [9]:
NQ;
Cl’
c*,
C3’
C4’
C5’
C6’
C7’
c,) = 0
K = (l_$)K
p=sp “V
m + @(S"VK
(2.18)
+ SwKcw)
+s&=+
(2.19)
with the following notation:
(2.10) c.
1
where the nonlinear operator is given by
g
h
= nonlinear =
coefficient,
gravitational
= enthalpy,
acceleration,
[LT-'1,
[L2Te2],
k kr
= permeability, [L2], = relative permeability,
P s T
= hydrodynamic pressure, = saturation, [LO], = temperature, ("K),
v
= velocity,
v
= specific volume,
K 4
= thermal conductivity, = porosity, IL'], = density, [ML-3 1, and
P 11 = viscosity,
(i = 1,...,8),
IL"], [ML-1T-2],
[LT-l], [L3Me1],
[ML-lTvl].
[MLT-3("K)-1],
1
118
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems
A numerical simulator of (2.10) can be written
(2.20)
where j denotes the nodal index. With proper initial condition, solution of (2.20) can be calculated provided that an accurate representation of the thermodynamics of steam-water systems is known. Such thermodynamics has been studied by Nguyen et al. [7] using nonequilibrium interpretation of the van der Waals equation of state which is referred to as the Riemann-Hugoniot In general, the new view of thermocatastrophe. dynamic evolution may be summarized as follows. We begin by identifying a location in thermodynamic space of Figure 6 designated by the numeral 1. This represents the initial condition or state at a physical point in our geothermal In this instance the point is in the simulator. hot water region. As the numerical simulation of this experiment evolves through time the pressure decreases dramatically and the thermodynamic state of the point of observation evolves along the trajectory defined by the curve (l-2). At the point denoted by the numeral 2, a second phase begins to appear. According to the catastrophe theory concept at this point, the specific volume of the steam phase is denoted by the thermodynamic state at point 4 and that of the water phase by the location 2. The trajectory continues along the line (2-3). During this interval of time the evolution of the specific volume of the water phase is described by the curve (2-3) and that of the steam phase by the curve (4-5). At the point 3, the thermodynamic trajectory describing the behavior at our obserSuddenly the vation point suddenly changes. state moves along the line (3-5) and only steam exists at this point in the system. Thereafter, the evolution of the state of the system at the observation point is described by the trajectory (5-6).
Figure 6.
While this description may appear quite straightIt proposes forward, it is certainly unorthodox. an evolution in water and steam densities within the two-phase region whereas current thinking assumes the steam and water densities do not change within the two phase region and are described by the two points 2 and 5 respectively as when the two fluids are maintained at equilibrium. The difference in our approach arises from our intention to model behavior in the twophase region when the fluids are not in equilibrium. Our approach also requires the existence of different pressures for water and steam withIn other words, the in the two-phase region. pi-essure of the steam associated with the point 4 is quite different than that of the water described by point 2 although these thermodynamic states would coexist in the proposed The existence of a higher thermodynamic model. pressure in the steam phase than the water phase is apparent to anyone who has watched a pan of Because the steam forms a bubble, boiling water. it must have a higher pressure than the surrounding water. The relationship between p, V van der Waals' equation:
F(p; v,T)
= (p +
3)(v
-
and T is the
b) - RT = 0
(2.21)
where a = 170.8 x lo3 N - m4/gi and b = 1.694 m3/gm for steam-water systems, R is the universal thermodynamic constant (8.314 J/OK-mole). The system of equations (2.20) and (2.21) has been solved by Nguyen and Pinder [8] with the finite-difference simulator shown in Figure 7.
Morphological Representation of SteamWater as Riemann-Hugoniot Catastrophe.
=%J r
Figure 7.
Flow Chart of the Geothermal Reservoir Simulator.
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems Details of the nonequilibrium scheme are described therein. Solution to the steam-water system for a standard set of initial and boundary conditions is plotted in Figures 8a and 8b. The results using the proposed simulator are more conservative than those obtained by equilibrium thermodynamics in place of the nonequilibrium scheme under the same boundary and initial conditions. Because the exact solution is not possible and experimental data are lacking, we cannot determine unequivocally which set of solutions is more A more subtle difference in the two accurate. approaches, equilibrium and nonequilibrium, involves the numerical treatment of the phase change. Because there is no sudden discontinuity in the derivatives of the nonlinear coefficient c.'s using the latter approach, the computational akpects of the problem are less troublesome, see 181 for a detailed analysis.
119
where 4 = porosity, [LO], = saturation of the a-phase,
[L"' Jr
sa a
= intrinsic mass density of the a-phase,
P
v
w31 I = pore velocity of the a-phase, [LT-l], and
a
Da = dispersion tensor, [L'T-l]. -0.
and e (pi) is an aggregate density of mass ex-
INITIAL PRESSURE ---__-_____________
change of species (i) into the u-phase from other phases. Equations (2.22) is subject to the constraints N M j sa=l (2.23) i=l Ci=l
ug= j
-
“ONEOUILIORIUU
T”ERYDDIWYl
(2.24)
Similar to equation (2.10), the generalized Darcy's law is used to represent va: Figure 8a.
906
Distribution of Nodal Pressure in Steam-Water System. -
NWEO”ILIBRI”Y
-
EO”lLlBRl”Y
k krcr " a =_=---.lJ,GSa (VP, - pagVD)
(2.25)
THERYOWNAYICS
ci = 1, . .. . M
THERYOOY)IIYICS
where k = absolute permeability z
of the solid matrix,
L21, kra= relative permeability uoi
= viscosity
of the a-phase,
of the a-phase,
= pressure of the a-phase
pcx D = depth, [L], and 0
.I2
.24
.?K
.46
60
g = gravitational
LE)(GT" , “ETERS Figure 8b. 2.4
Contrast in Enthalpy Profiles.
Multispecies/Multiphase voir [9]
By combining
[ML-~T-~],
[ML'~T-~],
acceleration,
-2 [LT 1.
(2.22) and (2.25) we obtain
Flow in Oil Reser-
Neglecting chemical reactions but allowing interphase mass transfer, the species balance equations governing the local average mass frac-
(2.26)
i=l
, .... N
H
where
tion tig of each species i = 1, . . .. N in every phase CL = l,..., M (excluding the immobile solid phase) are
ci
=
Cc=1
@ s,Pclo~ ,
i = 1,
. . .. N
(2.27)
120
V. V. Nguyen, E.F. Wood /Nonlinear
*.J
+SD
(2.29)
J _cY. = - (VP, - P'gVD)
(2.30)
I:
, i=
RT P=V-b-+wT$+i$) where a, b, y and 6 are
In practice the right-hand side of (2.26) is neglected provided that the distribution of each species (i) among phases is required to satisfy the local thermodynamic equilibrium conditions of Gibbs: f: = f: = ... = $
resources systems
For multispecies/multiphase flow in hydrocarbon reservoir, a practical equation of state known as the Developed Redlich-Kwong equation [lo] has recently been established, viz.
x = k kr,/uo mobility of the a-phase ..a
simulation of environmental
1, .... N
(2.31)
functions of temperature
and the distribution of species uy, i = l,...N. Let (2.37)
a2 = -
(RT+a(l+y))b P
(2.38)
where the functional forms for the fugacity coefficient fq of species in the a-phase are given, viz. (2.32)
a=1
,
..-,
a
CL =
= pacpa, OS, .. . . u;)
1,
(2.33)
__I6.a b3 P
then the normal form for (2.36) is the butterfly catastrophe:
F(v;
M ; i = 1, ... . N
Furthermore, to complete the system we also need M-l capillary pressure equations between (n,B) phases, and M equations of state P
a4 =
al,
a2,
a3,
a4) =
V5+ alV3
+aV+a 3
4
+ a2V3
=0
(2.41)
The geometric description of (2.41) is shown in Figures 9 and 10. The Developed Redlich-Kwong equation has been tested with field conditions for carbon dioxide flooding of oil reservoir. Depending on the position of (al, a2, a3, a,)
.. . . M
The system of relations (2.23) of M+l equations, (2.26) of N equations, (2.32) of N(M-1) equations, M-l capillary equations, (2.33) of M equations constitute the nonlinear operator N(s) for the state s = (WY, Sa, p,, fy, pa) of M(2N+3) If we define the bulk density of the unknowns. active phases to be M (2.34) P = 1 S/ a=1 then the specific volume of the system V = l/p. Given a thermodynamic pressure field p, neglecting the capillary pressure between phases, of all possible pa, it is important to determine the number M of stable phases existing at a temperature T(OK). If the state equation for p, p and V is known experimentally, then M and a P can be estimated from the number of real positive roots for V of F(v; p, T) = 0
(2.35)
O,IT.P)
0
and the principles of coexistence phases, see 171.
\
Figure 9.
A Geometric Description of the Butterfly Equation of State for Multispecies/Multiphase Hydrocarbon Systems.
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems
121
would like to thank Professor G.F. Pinder and Myron Allen for their comments on the multiphase flow problem in oil reservoir. REFERENCES
[ll
Thorn, R., Structural stability and morphogenesis, Trans. from the French by D.H. Fowler (Benjamin, Reading, Mass., 1975) 348 p.
[21 stewart, I., Applications of catastrophe theory to the physical sciences, Physica 20, (1981) 245-305. selected 131 Zeeman, c., Catastrophe theory: papers 1972-1977 (Addison-Wesley, Reading, Mass.) 675 p. 141 Nguyen, V.V. and Wood, E., On the Morphology of summer algae dynamics in non-stratified lakes, Ecological Modeling, 6 (1979) 117131. Figure 10.
Cross Sections of the Bifurcation Set of the Developed Redlich-Kwong Equation of State for Different Values of al(T,p) and a2(T,p). 1111
in the (p,T)-field, various combinations of coexisting phases can be realized 193. This information when used in parallel with the "umerical simulator of N(s) would facilitate the computational scheme in a manner similar to that of the geothermal problem. Results concerning this general class of multispecies/multiphase systems are intended for a future paper. CONCLUSIONS We have presented a general framework for the application of the elementary catastrophe theory to the simulation of environmental resources systems. Due to nonlinearity inherent in the dynamics of these systems, qualitative change in the structure of the state sometimes causes difficulties to the process of obtaining "umerThe normal forms of Thorn have ical solution. lent themselves to the investigation of the bifurcation set of external parameters which represents the predictable thresholds of the system and therefore help circumvent those diffiQuantitative results obtained for the culties. presented models have demonstrated many practical aspects of the proposed approach. It is hoped that this methodology may find itself useful in other engineering and applied sciences disciplines.
ISI Nguyen. V.V. and Wood, E.F., Bifurcation analysis of Rainfall infiltration, Water Resources Research, 17(l) (February 1981) 216-222. Nguyen, V.V. and Pinder, G.F., Is Geothermal simulation a catastrophe?, Proc. of the VI Annual Workshop on Geothermal Reservoir Engineering, Stanford University, California, (December 1980) 213-217. 171 Nguyen, V.V., Pinder, G.F. and J.F. Botha, Phenomenological interpretation of the thermodynamics of steam-water systems using catastrophe theory, J. Nonequilibrium Thermodynamics, 6 (1981), 285-295. T81 Nguyen, V.V. and Pinder, G.F., Geothermal reservoir simulation using nonequilibrium thermodynamics, Research Report 81-WR-2, Dept. of Civil Bngrg., Princeton University (July 1981), 35p (to appear in SPE J.). [91
Nguyen, V.V., Lectures on stability, catastrophes and chaos in simple dynamical systems, School of Appl. Sci: and Engrg., Princeton University, Fall and Spring 19801981, 430p. (monograph in preparation).
[lOI
Behar, E. and Jai", C., Application d'une equation d'btat 2 la simulation du comportment thermodynamique des fluides dans les gisement et dans les installations de surface, Rev. L'Institut Frangais du Petrole, 36(2) (1981), 173-181.
IllI
Brocker, T. and Lander, L., Differentiable germs and catastrophes, London Math. Sot. Lecture Notes 17, Cambridge Univ. Press, Cambridge, 1975, 179~.
ACKNOWLEDGEMENT This work was supported by the Department of Energy under Contract No. DE-ACO3-80F11489. The first part of the paper was presented at the IFIP Conference on Simulation and the Environment, Bangor, North Wales (September 19, 1979), when V.V. Nguyen was supported by the Wallace Fellowship of Princeton University. The authors
Mathematicsand Computersin SimulationXXIV
(1982) 114-121
North-Holland PublishingCompany
NONLINEAR
SIMULATION
OF ENVIRONMENTAL
RESOURCES SYSTEMS
V.V. NGUYEN and E.F. WOOD Water Resources
Program,
Department
of Civil Engineering,
Princeton
University,
Princeton,
NJ 08544,
U.S.A.
Nonlinearity, inherent in physical systems, can be simulated using elementary catastrophe theory. In this paper, we present a general framework for the application of this theory to environmental systems and demonstrate its practicality in a series of quantitative models consisting of the eutrophication of small lakes, rainfall infiltration process, geothermal simulation and the multispecies/multiphase flow in oil reservoir. INTRODUCTION The application of elementary catastrophe theory originated by Thorn [l] has evolved to a practical stage where the simulation of environmental systems can benefit from its unorthodox features. The state-of-the-art has been reviewed recently by Stewart [Z] and is well documented in [3].
Generally speaking, natural systems are often nonlinear and consequently exhibit qualitative changes in the structure of the state variables as related external or control parameters reach critical thresholds. Linear models, when employed to simulate these systems, fail to depict In order to propthis characteristic behavior. erly understand such systems, the geometric essence of the state must be brought forth according to conservation laws, constitutive axioms and intrinsic relationships between specific combination of external parameters and the state. When the behavior of the system is governed by the principle of structural stability, i.e., the description of the critical threshold remains invariant under perturbation of the governing energy field, the state-parameter relationships can be classified in terms of elementary catastrophes. For a discussion of these basic forms, we refer the reader to [l] "1. [3]. In this paper we explain how to employ elementary catastrophes to describe or numerically simulate environmental systems. A brief summary of a general framework of application is presented and is followed by a series of models. Two well-known forms are used in this presentation: the cusp and the butterfly catastrophes. Technical details of these models are given elsewhere, see 141, 151, [6], 171, and I&31. 1.
SUMMARY OF THE THEORETICAL FRAMEWORK _-
A mathematical description of a physical process subject to conservation laws and constitutive axioms usually has the form
0378-4754/82/0000-0000/$02.75
N(s;
a,b,c,...) = 0
(i.1)
where s is the state vector and ia,b,c,...} is the set of external or control parameters. The nonlinear operator N(e) consists of the time a (*): the local spatial changes, rate change, 5 such as a/a/ ~a dX
CL = (a 1'"" n lcll = c ui,
’ where x is a position vector,
o1 ), axa n
= ax
al a2 ax . ..ax 12
dt
,
and/or a combination of the latter.
When written as a discrete simulator, the form
*
cln
(1.1) has
c;. 2)
= [A]; + [B]
together with the initial conditions, s is the approximation of the state vector s at nodal ;oi;ts, [A] is the coefficient matrix containing a, ,c,... i and [I?]is the column vector of boundary and previous time level information. When the syste? is nonlinear, [A] will also be a function of s. In order to solve for sAin this situation, the relationship between s and [A] must be given a priori. For example, a; nodal points (31, the relationship between s i and [a,b,c,...) may be written as
P(sj; a,b,c,....) =
0
(1.3)
where P is a polynomial in the given variables. Notice that the time scale for relation (1.3) is assumed to be smaller than that of the evolution equation (1.2). Consequently, we can obtain dsj/dt from (;.3) (at least in principle and Since usually by somq iteration scheme). n s z (l.), (1.3) should be used in parallel to (1.2) in determining the approximate solution of (1.1). When we consider the principle of structural stability to establish (1.3), P(e) is usually best described by one of the seven ele-
0 1982 IMACS/North-Holland
V. V. Nguyen, E.F.
Wood /Nonlinear
simulation of environmental
mentary catastrophes of Thorn [I]. Each parameter a, b, c, . . .. may be regarded as a lumped parameter of other physical parameters of the dynamics under study. Let the normal form of P be
F(Sj; a,b,c,d,e) = 0
(1.4)
where the total number of lump parameters may not exceed five. From (1.4), the set of bifurcation values of {a,b,c,d,e> can be established by defining it to be the values that satisfy (1.5)
Essentially (1.5) denotes the coy?ditions by which the qualitative change in s at the smaller j time scale takes place with respect to change in As the external parameters cross m,d,e1. these thresholds in the control space, the method we use to solve (1.2) at the larger time scale must be adjusted accordingly. In other words, the loci of critical values of {a,b,c,d,e),i.e. the bifurcation set,partition the control space into different characteristic zones within which the corresponding behavior of the state s popsesses different structures. 2.
APPLICATIONS
We illustrate the above framework of utilizing the elementary catastrophe theory in the following analysis of environmental systems. 2.1
E; = F(x;
a,b)
(2.4)
By shortwhere E is a small positive constant. However, term, we mean two to three days periods. at a larger time scale, monthly or seasonally, the trend of growth in P is not steady, viz. dP/dt # 0
(2.5)
Therefore (2.4) may be used to specify the critical combination of chemical, mechanical and thermodynamic factors by which P experiences short-term blooms and crashes. The model of algae dynamics is presented with field data in Figures 1, 2, and 3. Chlorophyll a and ammonium nitrate are used as indicators of P and N ree' spectively.
300-
LAKE
555
SUWER 250
1972
-
3200. z
i
-SAMPLEDULUES
I u
no-
-sLorrcHANsE -II*py)owJGE
100 -
0
Let P denote the b&omass of phytoplankton representing the state s, by using the Monad and the Michaelis-Menten relations, we can show that
0
5
IO YAXIYUM
Figure 1. f2’
115
50.
Summer Algae Dynamics in Nonstratified Lakes 141
P@; fl’
resources systems
f3) = - (P3 + flP2 + f2P + f3) = 0
I5 Am
20
TEYPmmJRE.
25
30
-c
Phytoplankton Biomass as a Multivalued-Function of Air Temperature.
(2.1) where
fl
= fl
f2 = f2
(I,
Ne,
R)
(2.2.3)
(I,
Ne,
R)
(2.2b)
f3 = f3 (Ta, I) and
(2.2c)
I= light influx density, = chemical nutrient level in the lake e environment, R= endogeneous respiration factor, and Ta = air temperature. N
Using the transformation
x = P - f1/3, we can
reduce (2.1) to the cusp catastrophe form, i.e. F(x;
a,b) = - (x3 + a(N e,R)x + b(Ta,I))
3
(2.3)
According to the principle of structural stability [3],we may postulate the short-term dynamics of x by
Figure 2.
Control Space Trajectory for Lake 885.
116
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems where a, b, c are functions of oils.
When the
system conditions are known, i.e. 8,'s are deter1 mined a priori, one can find the critical rainfall rate r as a root of (2.9) at which the c upper layer of the porous medium assumes the state of ixmnediate ponding. This result implies different flow regimes associated with the wetting curve, see Figure 4. Denoting I to be the bulk infiltration volume, i.e., I = sdf where dz is the constant thickness of the moisture front, we obtain a morphological the variation of I in Figure 5.
descriotion of
-140 PLAINFIELD
Figure 3.
2.2
Morphological Representation of the state of P on a Cusp Catastrophe.
Physics of Rainfall Infiltration
[5]
The imbibition of water at the soil surface due to rainfall sometimes exhibits bifurcation behavior known as the Haines jumps. Using the physics of flow in the unsaturated zone, the combined equation of movement of air and water can be obtained according to the generalized Darcy's law. The integral form of this equation over a specific range of moisture content becomes (aHc + rt) 2
= P(s; ~1, ~2, p3, ~4)
REGION
where the state s is the normalized moisture content, H, is the effective capillary drive, r is the constant rainfall rate, t is the "shorttern?' time of precipitation, a is a gravitational characteristic constant, and p.'s are the para1 meters representing soil type, wetness condition, rainfall intensity and rainwater properties. Under commonly encountered conditions, the normal form of (2.6) may be approximated by 3
e3 + (0, +F)x
e5 + (e, +T)
(2.7)
are transformed representations t and (aHc + rt), respectively, and oils are
where x and E of
s
the transformed parameters of pi's_ the bifurcation 4(8
2
+2)3 r
From (1.5),
set is determined by
+ 27(8
4
+$2 r
= 0
(2.8)
which leads r
3
+ ar
2
+br+C=O
OF IMMEDIATE
?ONDING
(2.6) X.
dx Etx=X
SAND
(2.9)
Figure 4.
NDFWALIZED
MOISTURE CONTENT
Unsaturated Flow Regimes for Plainfield Sand.
V. V. Nguyen, E.F. Wood jNonlinear
117
simulation of environmental resources systems and
h,(e) = CAPILLARY DRIVE
’
a0
4%
I[
5
c2
c3
c4
=
“ah
ah (l-$)0 2 m ap
(l-$)P,
-4
2 -
4P
(2.12)
1 1 --c
5
-c
0 -c
6
=
7
k
1
K -$ +(; p,p,SvSJ(hv-h&e "V
k - +' ww
aT Kah (2.13)
'8
(2.14)
= $PV
subscripts m : solid phase, v : steam phase and w f liquid water phase,
h,(,&=SOIL
Figure 5.
h = t
PROPERTIES
Morphological Representation of Bulk Infiltration Volume over Capillary Hysteresis.
sv+s
(2.15)
(SvP,hv + SwP,hw)
w
(2.16)
=1
(2.17)
v = ; (S"U"V" + S"PWVJ 2.3
Geothermal Reservoir Simulation with Nonequilibrium Thermodynamics [7], [8]
The state s of a geothermal flow consists of In one dimension of pressure and enthalpy. space, the flow equations are represented by [9]:
NQ;
Cl’
c*,
C3’
C4’
C5’
C6’
C7’
c,) = 0
K = (l_$)K
p=sp “V
m + @(S"VK
(2.18)
+ SwKcw)
+s&=+
(2.19)
with the following notation:
(2.10) c.
1
where the nonlinear operator is given by
g
h
= nonlinear =
coefficient,
gravitational
= enthalpy,
acceleration,
[LT-'1,
[L2Te2],
k kr
= permeability, [L2], = relative permeability,
P s T
= hydrodynamic pressure, = saturation, [LO], = temperature, ("K),
v
= velocity,
v
= specific volume,
K 4
= thermal conductivity, = porosity, IL'], = density, [ML-3 1, and
P 11 = viscosity,
(i = 1,...,8),
IL"], [ML-1T-2],
[LT-l], [L3Me1],
[ML-lTvl].
[MLT-3("K)-1],
1
118
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems
A numerical simulator of (2.10) can be written
(2.20)
where j denotes the nodal index. With proper initial condition, solution of (2.20) can be calculated provided that an accurate representation of the thermodynamics of steam-water systems is known. Such thermodynamics has been studied by Nguyen et al. [7] using nonequilibrium interpretation of the van der Waals equation of state which is referred to as the Riemann-Hugoniot In general, the new view of thermocatastrophe. dynamic evolution may be summarized as follows. We begin by identifying a location in thermodynamic space of Figure 6 designated by the numeral 1. This represents the initial condition or state at a physical point in our geothermal In this instance the point is in the simulator. hot water region. As the numerical simulation of this experiment evolves through time the pressure decreases dramatically and the thermodynamic state of the point of observation evolves along the trajectory defined by the curve (l-2). At the point denoted by the numeral 2, a second phase begins to appear. According to the catastrophe theory concept at this point, the specific volume of the steam phase is denoted by the thermodynamic state at point 4 and that of the water phase by the location 2. The trajectory continues along the line (2-3). During this interval of time the evolution of the specific volume of the water phase is described by the curve (2-3) and that of the steam phase by the curve (4-5). At the point 3, the thermodynamic trajectory describing the behavior at our obserSuddenly the vation point suddenly changes. state moves along the line (3-5) and only steam exists at this point in the system. Thereafter, the evolution of the state of the system at the observation point is described by the trajectory (5-6).
Figure 6.
While this description may appear quite straightIt proposes forward, it is certainly unorthodox. an evolution in water and steam densities within the two-phase region whereas current thinking assumes the steam and water densities do not change within the two phase region and are described by the two points 2 and 5 respectively as when the two fluids are maintained at equilibrium. The difference in our approach arises from our intention to model behavior in the twophase region when the fluids are not in equilibrium. Our approach also requires the existence of different pressures for water and steam withIn other words, the in the two-phase region. pi-essure of the steam associated with the point 4 is quite different than that of the water described by point 2 although these thermodynamic states would coexist in the proposed The existence of a higher thermodynamic model. pressure in the steam phase than the water phase is apparent to anyone who has watched a pan of Because the steam forms a bubble, boiling water. it must have a higher pressure than the surrounding water. The relationship between p, V van der Waals' equation:
F(p; v,T)
= (p +
3)(v
-
and T is the
b) - RT = 0
(2.21)
where a = 170.8 x lo3 N - m4/gi and b = 1.694 m3/gm for steam-water systems, R is the universal thermodynamic constant (8.314 J/OK-mole). The system of equations (2.20) and (2.21) has been solved by Nguyen and Pinder [8] with the finite-difference simulator shown in Figure 7.
Morphological Representation of SteamWater as Riemann-Hugoniot Catastrophe.
=%J r
Figure 7.
Flow Chart of the Geothermal Reservoir Simulator.
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems Details of the nonequilibrium scheme are described therein. Solution to the steam-water system for a standard set of initial and boundary conditions is plotted in Figures 8a and 8b. The results using the proposed simulator are more conservative than those obtained by equilibrium thermodynamics in place of the nonequilibrium scheme under the same boundary and initial conditions. Because the exact solution is not possible and experimental data are lacking, we cannot determine unequivocally which set of solutions is more A more subtle difference in the two accurate. approaches, equilibrium and nonequilibrium, involves the numerical treatment of the phase change. Because there is no sudden discontinuity in the derivatives of the nonlinear coefficient c.'s using the latter approach, the computational akpects of the problem are less troublesome, see 181 for a detailed analysis.
119
where 4 = porosity, [LO], = saturation of the a-phase,
[L"' Jr
sa a
= intrinsic mass density of the a-phase,
P
v
w31 I = pore velocity of the a-phase, [LT-l], and
a
Da = dispersion tensor, [L'T-l]. -0.
and e (pi) is an aggregate density of mass ex-
INITIAL PRESSURE ---__-_____________
change of species (i) into the u-phase from other phases. Equations (2.22) is subject to the constraints N M j sa=l (2.23) i=l Ci=l
ug= j
-
“ONEOUILIORIUU
T”ERYDDIWYl
(2.24)
Similar to equation (2.10), the generalized Darcy's law is used to represent va: Figure 8a.
906
Distribution of Nodal Pressure in Steam-Water System. -
NWEO”ILIBRI”Y
-
EO”lLlBRl”Y
k krcr " a =_=---.lJ,GSa (VP, - pagVD)
(2.25)
THERYOWNAYICS
ci = 1, . .. . M
THERYOOY)IIYICS
where k = absolute permeability z
of the solid matrix,
L21, kra= relative permeability uoi
= viscosity
of the a-phase,
of the a-phase,
= pressure of the a-phase
pcx D = depth, [L], and 0
.I2
.24
.?K
.46
60
g = gravitational
LE)(GT" , “ETERS Figure 8b. 2.4
Contrast in Enthalpy Profiles.
Multispecies/Multiphase voir [9]
By combining
[ML-~T-~],
[ML'~T-~],
acceleration,
-2 [LT 1.
(2.22) and (2.25) we obtain
Flow in Oil Reser-
Neglecting chemical reactions but allowing interphase mass transfer, the species balance equations governing the local average mass frac-
(2.26)
i=l
, .... N
H
where
tion tig of each species i = 1, . . .. N in every phase CL = l,..., M (excluding the immobile solid phase) are
ci
=
Cc=1
@ s,Pclo~ ,
i = 1,
. . .. N
(2.27)
120
V. V. Nguyen, E.F. Wood /Nonlinear
*.J
+SD
(2.29)
J _cY. = - (VP, - P'gVD)
(2.30)
I:
, i=
RT P=V-b-+wT$+i$) where a, b, y and 6 are
In practice the right-hand side of (2.26) is neglected provided that the distribution of each species (i) among phases is required to satisfy the local thermodynamic equilibrium conditions of Gibbs: f: = f: = ... = $
resources systems
For multispecies/multiphase flow in hydrocarbon reservoir, a practical equation of state known as the Developed Redlich-Kwong equation [lo] has recently been established, viz.
x = k kr,/uo mobility of the a-phase ..a
simulation of environmental
1, .... N
(2.31)
functions of temperature
and the distribution of species uy, i = l,...N. Let (2.37)
a2 = -
(RT+a(l+y))b P
(2.38)
where the functional forms for the fugacity coefficient fq of species in the a-phase are given, viz. (2.32)
a=1
,
..-,
a
CL =
= pacpa, OS, .. . . u;)
1,
(2.33)
__I6.a b3 P
then the normal form for (2.36) is the butterfly catastrophe:
F(v;
M ; i = 1, ... . N
Furthermore, to complete the system we also need M-l capillary pressure equations between (n,B) phases, and M equations of state P
a4 =
al,
a2,
a3,
a4) =
V5+ alV3
+aV+a 3
4
+ a2V3
=0
(2.41)
The geometric description of (2.41) is shown in Figures 9 and 10. The Developed Redlich-Kwong equation has been tested with field conditions for carbon dioxide flooding of oil reservoir. Depending on the position of (al, a2, a3, a,)
.. . . M
The system of relations (2.23) of M+l equations, (2.26) of N equations, (2.32) of N(M-1) equations, M-l capillary equations, (2.33) of M equations constitute the nonlinear operator N(s) for the state s = (WY, Sa, p,, fy, pa) of M(2N+3) If we define the bulk density of the unknowns. active phases to be M (2.34) P = 1 S/ a=1 then the specific volume of the system V = l/p. Given a thermodynamic pressure field p, neglecting the capillary pressure between phases, of all possible pa, it is important to determine the number M of stable phases existing at a temperature T(OK). If the state equation for p, p and V is known experimentally, then M and a P can be estimated from the number of real positive roots for V of F(v; p, T) = 0
(2.35)
O,IT.P)
0
and the principles of coexistence phases, see 171.
\
Figure 9.
A Geometric Description of the Butterfly Equation of State for Multispecies/Multiphase Hydrocarbon Systems.
V. V. Nguyen, E.F. Wood /Nonlinear simulation of environmental resources systems
121
would like to thank Professor G.F. Pinder and Myron Allen for their comments on the multiphase flow problem in oil reservoir. REFERENCES
[ll
Thorn, R., Structural stability and morphogenesis, Trans. from the French by D.H. Fowler (Benjamin, Reading, Mass., 1975) 348 p.
[21 stewart, I., Applications of catastrophe theory to the physical sciences, Physica 20, (1981) 245-305. selected 131 Zeeman, c., Catastrophe theory: papers 1972-1977 (Addison-Wesley, Reading, Mass.) 675 p. 141 Nguyen, V.V. and Wood, E., On the Morphology of summer algae dynamics in non-stratified lakes, Ecological Modeling, 6 (1979) 117131. Figure 10.
Cross Sections of the Bifurcation Set of the Developed Redlich-Kwong Equation of State for Different Values of al(T,p) and a2(T,p). 1111
in the (p,T)-field, various combinations of coexisting phases can be realized 193. This information when used in parallel with the "umerical simulator of N(s) would facilitate the computational scheme in a manner similar to that of the geothermal problem. Results concerning this general class of multispecies/multiphase systems are intended for a future paper. CONCLUSIONS We have presented a general framework for the application of the elementary catastrophe theory to the simulation of environmental resources systems. Due to nonlinearity inherent in the dynamics of these systems, qualitative change in the structure of the state sometimes causes difficulties to the process of obtaining "umerThe normal forms of Thorn have ical solution. lent themselves to the investigation of the bifurcation set of external parameters which represents the predictable thresholds of the system and therefore help circumvent those diffiQuantitative results obtained for the culties. presented models have demonstrated many practical aspects of the proposed approach. It is hoped that this methodology may find itself useful in other engineering and applied sciences disciplines.
ISI Nguyen. V.V. and Wood, E.F., Bifurcation analysis of Rainfall infiltration, Water Resources Research, 17(l) (February 1981) 216-222. Nguyen, V.V. and Pinder, G.F., Is Geothermal simulation a catastrophe?, Proc. of the VI Annual Workshop on Geothermal Reservoir Engineering, Stanford University, California, (December 1980) 213-217. 171 Nguyen, V.V., Pinder, G.F. and J.F. Botha, Phenomenological interpretation of the thermodynamics of steam-water systems using catastrophe theory, J. Nonequilibrium Thermodynamics, 6 (1981), 285-295. T81 Nguyen, V.V. and Pinder, G.F., Geothermal reservoir simulation using nonequilibrium thermodynamics, Research Report 81-WR-2, Dept. of Civil Bngrg., Princeton University (July 1981), 35p (to appear in SPE J.). [91
Nguyen, V.V., Lectures on stability, catastrophes and chaos in simple dynamical systems, School of Appl. Sci: and Engrg., Princeton University, Fall and Spring 19801981, 430p. (monograph in preparation).
[lOI
Behar, E. and Jai", C., Application d'une equation d'btat 2 la simulation du comportment thermodynamique des fluides dans les gisement et dans les installations de surface, Rev. L'Institut Frangais du Petrole, 36(2) (1981), 173-181.
IllI
Brocker, T. and Lander, L., Differentiable germs and catastrophes, London Math. Sot. Lecture Notes 17, Cambridge Univ. Press, Cambridge, 1975, 179~.
ACKNOWLEDGEMENT This work was supported by the Department of Energy under Contract No. DE-ACO3-80F11489. The first part of the paper was presented at the IFIP Conference on Simulation and the Environment, Bangor, North Wales (September 19, 1979), when V.V. Nguyen was supported by the Wallace Fellowship of Princeton University. The authors