Dirichlet conditions are required on boundaries where the hydraulic head on the boundary doesn’t depend on the ow conditions in the aquifer. It will generally be where the aquifer is in direct contact with free water, such as a river, a lake or a sea. Along this contact between the aquifer and the river(or lake, sea. . . ), the hydraulic head is constant and prescribed by the water elevation in the river. The river can either feed or drain the aquifer. Of course, the water level in the river can change along its course, but the river still prescribes the hydraulic head along the boundary.
It is pretty obvious to imagine how we can take into account this type of boundary condition in the kriging system : by discretizing the continuous boundary into a nite number of data points which will be assigned the prescribed head value(s). We will simply add these points to the data points provided by water table measurements.
This is the Neumann condition in hydrogeology. According to the Darcy law, prescribing the head gradient normal to the boundary, @h @n , is indeed the same as prescribing the ux −T @h @n on this boundary, provided T is known. There are two distinct conditions of prescribed ux :
• The no ow boundaries : @h @n = 0. For example, in a 2D ow, the contact between an aquifer and neighboring impervious formations.
• The prescribed ux with a value dierent than 0. For example, runo water entering an aquifer along a boundary.
To implement this new type of kriging, I used the GSLIB open source code, developed in Stanford University, and documented in Deutsch and Journel (1998). This code is popular both among researchers and professionals, and it is particularly used in Waterloo Hydrogeologic softwares. The code implementation can be divided in several units. The algorithms won’t be detailed in this section, but they are provided in the appendixes. This section starts with a common kriging issue before presenting the main algorithms introduced.
The kriging theory is always derived as if all the N data points were used in the estimation ; this is the so-called global neighborhood case. In practice, N may be too large to allow computation and a moving neighborhood has to be used, including only a subset of the data for the estimation of each grid node (see gure 2.2). Formally, this does not change anything for a grid node taken in isolation : the content of the sampled set of points is just dierent. However, it may alter the relationships between estimates at dierent grid nodes and introduce spurious discontinuities.
Adding the Boundary Conditions Data
There are two input les in GSLIB : the rst one describes all the parameters the main algorithm needs and the second one provides the coordinates of the data points and the values of the studied regionalized variable, the hydraulic head in our case. The method chosen to introduce the boundary conditions data is to discretize the boundary lines into points. The user will input the coordinates and the head value for each point that is a boundary segment end. However, the program has to know if a point is a boundary segment end or not. And in the former case, it also has to know if it is a prescribed head or a prescribed ux condition. In order to solve this issue, a new input parameter has been created, named kod. This kod is set :
• to 0 for the measured data points.
• to a positive integer value for a prescribed head segment end.
• to a negative integer value for a prescribed ux segment end.
Points that are both ends of the same segment have the same kod. More generally, if a boundary is represented by a broken line, a point will be placed at each direction shift, and all these shift points should have the same kod, as they are part of the same boundary. Otherwise, two points representing two dierent boundary lines should have a dierent kod.
With these added points, we now have all the needed data for the Kriging under Boundary Conditions process, but we still don’t have really our discretized boundaries. In order to provide these, a simple algorithm will read the input data, detect the boundary points (the ones whose kod is not 0) and create several points between the two ends of the segment by linearly interpolating the coordinates (and the head value for a prescribed head segment). We have to notice here that, for this subroutine to work properly, the two ends of a boundary segment, or two consecutive shift points of a broken line, have to follow one another in the list of input data points. The data set now consists of both the real data points and the points representing the discretized boundary lines. Finally, another algorithm scans the data set again, selects the prescribed ux points, and adds the dummy points in order for them to make a segment perpendicular to the boundary line they are representing. For a dened boundary, this is done using the previous and the following boundary points to compute the local slope of the boundary and setting the slope of the dummy points segment to make it perpendicular to the boundary (cf. Fig. 2.1). For the segment ends, the local slope is computed using the end point itself and its closest neighbor in the segment.
Table of contents :
List of Figures
1.1.1 Random Variable
1.1.2 Random Function
1.1.3 Other Denitions
1.1.4 Stationarity Hypothesis
1.1.5 A Useful Result
1.2 Structural Analysis
1.2.1 Experimental Variogram
1.2.2 Variogram Characteristics
1.2.3 Variogram Examples
1.3 Kriging Basics
1.4 Universal Kriging or Kriging with a Trend Model
1.4.1 Drift Terms
1.4.2 Linearity Constraint
1.4.3 Authorization Constraint
1.4.4 Unbiasedness Constraints
1.4.5 Optimality Constraint
1.4.6 Solving the Kriging Equations
1.5 Multivariate Geostatistics : Cokriging
2 Kriging under Boundary Conditions
2.1 Boundary Conditions in Hydrogeology
2.1.1 Prescribed Head
2.1.2 Prescribed Flux
2.2 The Kriging under Boundary Conditions System
2.2.1 Linearity Constraint
2.2.2 Authorization Constraint
2.2.3 Unbiasedness Constraints
2.2.4 Optimality Constraint
2.3 Code Implementation
2.3.1 Kriging Neighborhood
2.3.2 Adding the Boundary Conditions Data
2.3.3 Discretization Parameters
2.3.4 Singularity Conditions
2.3.5 Constant Flux Conditions
2.3.6 Cubic Variogram
2.3.7 Filling the Kriging under Boundary Conditions Matrix
2.4 Application of Kriging under Boundary Conditions
2.4.1 Comparison between Universal Kriging and Kriging under Boundary
2.4.2 Undened boundaries
2.4.3 Constant Flux Boundaries
2.4.4 River and Inside Constant Flux
2.4.5 Screen Eect
3 (Co-)Kriging with Multivariate Structural Analysis
3.1 Further Geostatistical Denitions
3.1.1 Intrinsic Random Function of order k (IRF-k)
3.1.2 Generalized Covariance
3.2 Kriging the Head using Transmissivity Knowledge
3.2.1 Hydrogeologic Context
3.2.2 The Stochastic Equation Z = Y
3.2.3 Covariance Model
3.2.4 Covariance Choice and Code Implementation .
3.3 Cokriging Head and Log Transmissivity
3.3.2 The Cokriging System
3.3.3 The Cokriging under Boundary Conditions System
3.4 Inverse Problem
3.4.1 The Inverse Problem Kriging System
3.4.2 The Bias in Lognormal Kriging
3.5.1 Code Implementation