BME studies of stochastic differential equations representing physical laws-part II

BME STUDIES OF STOCHASTIC DIFFERENTIAL EQUATIONS REPRESENTING PHYSICAL LAWS -PART IIM.L. Serre and G. ChristakosEnvironmental Modelling Program, Department of Environmental Sciences & Engineering School of Public Health, University of North Carolina, Chapel Hill, NC 27599-7400, USA Classical Geostatistics methods have been designed to use mainly statistical knowledge about natural variables and they lack the ability to incorporate important forms of knowledge like physical laws and scientific theories into the mapping process. On the other hand, the powerful and versatile Bayesian maximum entropy (BME) method of Modern Geostatistics can accomplish such a task, rigorously and efficiently. In this work, BME is used to incorporate the Darcy law of subsurface hydrology in the spatial mapping of a hydraulic head field. The hydraulic map thus obtained is physically meaningful as well as numerically more accurate than that obtained using classical methods (e.g., kriging). Moreover, taking advantage of the Darcy law the BME hydraulic head mapping can involve other related soil properties, like hydraulic conductivity. The approach leads to very accurate hydraulic head solutions, and may be also applied to study the inverse problem, in which one seeks to estimate hydraulic conductivity from hydraulic head measurements.1. INTRODUCTIONThe Bayesian Maximum Entropy (BME) method of Modern Geostatistics is a method which offers a rigorous framework to account for a large class of general and specificatory knowledge bases [1-6]. The general knowledge G includes physical laws, scientific theories, empirical relationships, and statistical moments. The specificatory knowledge S includes hard data (exact measurement) and soft data (such as intervals, probability assessments, uncertain observations and fuzzy sets) that are specific to the mapping situation. Previous studies have shown that the BME method is very useful and numerically efficient in a number of practical situations in which various forms of soft data are available [5, 6]. In the first part of this work ([4]) we discussed two methods used by BME analysis in order to account for physical knowledge in the form of physical laws. In this work, the BME formalism is used to incorporate the Darcy law of groundwater flow in porous media, which relates changes in hydraulic head with soil properties such as the hydraulic conductivity or resistivity. Darcy law is considered part of the general knowledge G which leads to a prior G-based pdf that describes the joint spatial distribution of hydraulic head and resistivity. Then, the specificatory knowledge S is used to update the prior pdf leading to the posterior pdf of the hydraulic head map. Numerical examples show that, by accounting for the Darcy law BME provides predictions that are physically meaningful and numerically more accurate than classical data analysis that does not account for the physical law.2. INCORPORATING DARCY LAW IN BME MAPPINGConsider the Darcy law of one-dimensional flow usually expressed as K(s)∇H(s)=−q(s), where the specific discharge q(s) is considered deterministic while the hydraulic head H(s) and the hydraulic conductivity K(s) are random fields. It is convenient to define the hydraulic resistivity R(s)=1/K(s) and rewrite Darcy law asdH(s)/ds=−q(s)R(s).(1)Let H0 be the known value of H(s) at the spatial origin s=0, let Hmap=[H1,...,H nH]T be thevector of random variables representing the hydraulic head H(s) at points s1,...,sn H, and assumethat R data =[R n H +1,...,R n ]T is the vector of random variables representing the hydraulic resistivity R (s ) at points s n H +1,...,s n . At the prior stage of the BME analysis we seek to obtain the prior pdff G (H map ,R data ;s 1...s n ) given general knowledge G (which includes the Darcy law),where H map and R data represent a realization of the random vectors H map and R data ,respectively. The general knowledge G is transformed in a suitable set of moment equations as followsd H map ∫d R data ∫g α(H map ,R data )f G (H map ,R data ;s 1...s n )=g α (α=1,...,N α), where g α isa set of functions chosen such that their stochastic expectations g α are either known statistical moments or they can be inferred from Darcy law. For example, we assume that the hydraulic resistivity has a known mean m R (s )=R (s ), non-centered covariance C R (s ,s ')=R (s )R (s '), and 4th-order non-centered moment m R ,4(s )=R 4(s ). These moments are incorporated in the general knowledge by letting g i (R i ;s i )=R i and g i (s i )=m R (s i ), i ∈J R ={n H +1,...,n };g ij (R i ,R j ;s i ,s j )=R i R j and g ij (s i ,s j )=C R (s i ,s j ) for i ,j ∈J R ; and finally g n +i (R i ;s i )=R i 4and g n +i (s i )=m R ,4(s i ) for n +i ∈J R . On the other hand the moments describing the hydraulic head are not known a priori, but they may be inferred from Darcy law, Eq. (1).Using the Darcy law we may obtain expressions for the mean m H (s )=H (s ) and non-centered covariance function C H (s ,s ')=H (s )H (s ') of the hydraulic head, as well as the non-centered cross-covariance function C RH (s ,s ')=R (s )H (s '). Clearly, if we take the stochastic expectation of Eq. (1) and solve with respect to the boundary condition H 0 we obtain m H (s )=H o −du 0s∫q (u )m R (u ), which is readily incorporated into G by choosingg i (H i ;s i )=H i and g i (s i )=m H (s i ) for i ∈J H , where J R ={n H +1,...,n }. Similarly,multiplying Eq. (1) by R (s '), taking the stochastic expectation and solving, we obtain C HR (s ,s ')=H 0m R (s ')−duq (u )C R (u ,s ')0s∫, which we incorporate in G by choosingg ij (H i ,R j ;s i ,s j )=H i R j and g ij (s i ,s j )=C HR (s i ,s j ) for i ∈J H , j ∈J R . Other moments may be obtained in a similar way. For example, multiplying Eq. (1) by H (s ') and taking the stochasticexpectation leads C H (s ,s ')=du 0s∫du '0s '∫q (u )q (u ')C R (u ,u ')−H 02+H 0[m H (s )+m H (s ')], whichwe also incorporate in G by choosing g ij (H i ,H j ;s i ,s j )=H i H j and g ij (s i ,s j )=C H (s i ,s j ) for i ,j ∈J H . The selected g α functions, summarized in Table 1, represent a set of constraints which incorporates the general knowledge consisting of statistical moments and the Darcy ing these constraints under the epistemic goal of information maximization leads to the following solution for the prior pdf [3]f G (H map ,R data ;s 1...s n )=Z −1exp[µi i ∈J H ∑H i +µi i ∈J R ∑R i +µij j ∈J H∑i ∈J H ∑H i H j +µij j ∈J R∑i ∈J R ∑R i R j+µij j ∈J R∑i ∈J H ∑H i R i +µi +n i ∈J R∑R i 4](2)where the Lagrange coefficients µi and µij are calculated by solving the system of equations obtained from Table (1).In the posterior stage of the BME method, the specificatory knowledge S is used to update the prior pdf and lead to the posterior pdf f K (H k ) which completely describes the hydraulicTable 1: General knowledge constraintsαg αα=i , i ∈J g (H ;s )=H g (s )=m (s )α=i , i ∈J Rg i (R i ;s i )=R i g i (s i )=m R (s i )α=(ij ), i ∈J H , j ∈J Hg ij (H i ,H j ;s i ,s j )=H i H j g ij (s i ,s j )=C H (s i ,s j )α=(ij ), i ∈J R , j ∈J Rg ij (R i ,R j ;s i ,s j )=R i R j ;g ij (s i ,s j )=C R (s i ,s j )α=(ij ), i ∈J H , j ∈J Rg ij (H i ,R j ;s i ,s j )=H i R j ;g ij (s i ,s j )=C HR (s i ,s j )α=n +i , i ∈J Rg n +i (R i ;s i )=R i 4;g n +i (s i )=m R ,4(s i )head distribution. Let H hard T =[H 1,...,H m Hh ], H soft T=[H m Hh +1,...,H m H ], and H k represent H (s )at the hard data points, soft data points and estimation points, respectively, such that H map T =[H hard T H soft T H k ]. Similarly let R hard T =[R 1,...,R m Rh ] and R soft T =[R m Rh +1,...,R m R ]represent R (s ) at the hard and soft data points, respectively, such that R data T =[R hard T R soft T ]. For soft data of the interval type, the BME posterior pdf is given by ([3,4])f K (H k )=A −1d H soft l H u H ∫d R soft l Ru R∫f G (H map ,R data ), where l H and u H are the lower and upperbounds of the soft data for the hydraulic head such that l H ≤H soft ≤u H , and l R and u R are the lower and upper bounds of the soft data for the hydraulic resistivity so that l R ≤R soft ≤u R .3. SIMULATED CASE STUDYIn this case study we consider the R (s ) profile shown in Fig 1a. These values were generated at grid-points separated by a distance of 0.5 m , using the exponential covariance model c R (s ,s ')=c R ,0exp(−|s '−s |/a R ) with c R ,0=1 (sec/mm )2 and a R =6 m , and a mean of m R =2 sec/mm . The histogram of R (s ) values (Fig. 1b) is approximately symmetric with akurtosis coefficient of m R ,4/c R ,02≅3, while the corresponding K (s )=1/R (s ) values have an asymmetric distribution (Fig. 1c). The hydraulic head H (s ) values calculated by solving Darcy law, Eq. (1), with H 0=100 m and q (s )=q =0.05 mm /sec are shown in Fig. 2a. The mean hydraulic head m H (s )=H o −q m R s is shown with a dashed line in Fig. 2a, and the histogram of hydraulic head fluctuation H (s )−m H (s ) is shown with plain lines (Fig. 2b) at points corresponding to s =5 m , s =25 m and s =75 m . These histograms were obtained by using 2,000 H (s ) realizations, and they show that the variance of H (s ) increases as s increases. In Fig. 2b (dotted line) the prior pdf's for the hydraulic head at points s =5 m , s =25 m and s =75 m . Note that these are in excellent agreement with the experimental histogram calculated from the data.In a first test case, we consider the estimation of H (s )−m H (s ) using sparse measurements for H (s ) and R (s ). The selected hard and soft (interval) data are shown with triangles and error bars, respectively, in Fig 1a for R (s ), and in Fig 2a for H (s ). In Fig 3a the estimated head fluctuation profile was derived from simple kriging (SK) using only hard head data; the actual head fluctuation profile is also shown for comparison. Note the poor SK estimates at unobserved locations. The head fluctuation profile in Fig. 3b was obtained from BME using hard and soft (interval) head data as well as the Darcy law. Note that the BME estimated profile5010015020000.511.522.533.544.5Distance s (in meter)R (s ) (i n s e c /m m )(a)K (mm/sec)F r e q u e n c y(c)R (sec/mm)F r e q u e n c y(b)Figure 1: (a) Simulated profile for the hydraulic resistivity R (s ). Histograms of (b) the hydraulicresistivity R (s ), and (c) the hydraulic conductivity K (s )=1/R (s ).respectively. (b) H (s ) histogram at different spatial locations s .is a substantial improvement over the classical SK method. Indeed, by being able to incorporate the Darcy law, the BME method allowed to account for hard and soft (interval) hydraulic conductivity data, as well.In a second case study we consider the problem of estimating the hydraulic resistivity R (s )using sparse measurements of R (s ) and H (s ). This is sometimes called the inverse problem.While it is possible to obtain accurate measurements of the hydraulic head H (s ), it is difficult to directly measure the hydraulic resistivity. Therefore, in case study 2, we consider that the measurements for H (s ) are the hard data shown with triangles in Fig. 4a, while the information for R (s ) consists of the soft interval data shown with error bars in Fig. 4b. For illustration purposes we calculate BME estimates for 3 different mapping situations. In the first situation we assume that the specificatory knowledge consists of only the soft data for R (s ), and we show the corresponding BME profile for R (s ) in Fig. 4b. Note that this type of soft information commonly arises when estimating hydraulic resistivities, and the BME method provides a rigorous framework to account for such soft information. In the second mapping situation we assume that the specificatory knowledge includes both the soft data for R (s ), as well as the hard data for H (s ) shown with triangles in Fig. 4a. Using this information we obtain the profile shown in Fig. 4c, which clearly provides a substantial improvement over that of Fig. 4b. Note that the profile of Fig. 4c respects the soft data but also reproduces the important features of theH y d r a u l i c H e a d F l u c t u a t i o n (i n m e t e r )(triangles); (b) BME using hard and soft (interval) head data and the Darcy law. Dotted linesrepresent the actual head fluctuation profile.resistivity R (s ) obtained from: (b) BME using R (s ) data; (c) BME using data for both R (s ) and H (s ); and (d) BME using only H (s ) data. The dotted line represents the true R (s ) profile.true resistivity profile shown in dotted line. This unique feature of the BME approach is due to the fact that BME accounts for the Darcy law which relates H (s ) and R (s ), leading to substantial improvements in the estimation of the hydraulic resistivity. Finally, in Fig. 4d we show the R (s ) profile using only hard H (s ) data (triangles in Fig. 4a). Note that in this case the important features of the hydraulic resistivity profiles are well preserved.ACKNOWLEDGMENTS: This work was supported by grants from the National Institute of Environmental Health Sciences (Grant no. P42 ES05948-02) the Department of Energy (Grant no. DE-FC09-93SR18262), and the Army Research Office (Grant no. DAAG55-98-1-0289), and by the Computational Science Graduate Fellowship Program of the Office of Scientific Computing.REFERENCES1.CHRISTAKOS, G., "A Bayesian/maximum-entropy view to the spatial estimationproblem". Mathematical Geology, 22(7), 763-776 (1990).2.CHRISTAKOS, G.- Random Field Models in Earth Sciences. Academic Press,San Diego, CA, 1992.3.CHRISTAKOS, G. and X. LI - "Bayesian maximum entropy analysis and mapping:farewell to kriging estimators?" Mathematical Geology,30(4), 435-462 (1998).4.CHRISTAKOS, G, D.T. HRISTOPULOS and SERRE, M.L. - BME studies of stochasticdifferential equation representing physical laws -part I". In Proc. of the 5th Annual Confer.of the Internat. Assoc.for Mathematical Geology, Trodheim, Norway (1999).5.SERRE, M. L., P. BOGAERT and G. CHRISTAKOS. "Latest Computational Results inSpatiotemporal Prediction Using the Bayesian Maximum Entropy Method". In Proc. of the 4th Annual Confer. of the Internat. Assoc. for Mathematical Geology, Eds. Buccianti A., Nardi G., Potenza R., De Frede Editore, Napoli, 117-122 (1998).6.SERRE, M. L. and CHRISTAKOS, G.- "Modern geostatistics: computational BMEanalysis in the light of uncertain physical knowledge - the Equus Beds study". Journal of Stochastic Environmental Research and Risk Analysis, 13(1), 1-26 (1999)。