Since a geostatistical analysis was
used to estimate the complete potentiometric surface, a semi-variogram
was constructed from the sample data.
Investigation in various directions revealed that a strong trend was
present on the surface. This is not surprising with this type of variable and
indicates a strong flow direction, generally to the north-east
within the study area. Removing the
trend with a common Polynomial Trend Surface technique, allows semi-variograms
to be constructed for the residuals from the trend
surface. These are
seen to be isotropic, allowing a semi-variogram model to be fitted to
the experimental semi-variogram. The
model used is a Spherical semi-variogram with a range of influence of 50 miles,
a sill on the spherical component of 46,000 (ft)@ and a nugget effect of 11,000 (ft)@. This model was confirmed by a
cross validation study using Universal Kriging to include the form of the known
trend over the study area. A fuller
description of this study may be found in Harper et
al., 1986.
Once a semi-variogram model has been fitted
to the sample data over the study area, this model - combined with the original sample data - may be used to estimate values at selected locations
within that area. To produce a map of the.area it is usual to estimate values
at the nodes of a square grid laid over the area of interest. The size of the grid is usually a compromise
between the amount of computer time taken for the estimation and the amount of detail which is necessary on the final map. Figure 2 shows a contour map which was
constructed on a 5 mile grid laid across the study
area indicated in Figure 1. The map also shows the approximate site of the
planned high level nuclear waste repository. The solid line which has
been superimposed on this map is the 'most likely" travel path of a
freed radio-nuclide if the repository was breached. That is, this line follows the steepest slope
of
the potentiometric surface away from the site location. The path comes to a stop in an estimated
local low point in the Amarillo basin. This would generally be quoted as the travel path in a
report. There are two problems with this
prediction which we shall discuss here.
Uncertainty
on the estimated surface
When a value is estimated by a
geostatistical method, like Universal Kriging, two figures are produced. A Kriging technique produces an estimator
formed from a linear combination of the surrounding sample values. The particular linear combination used is
that which minimises the "estimation variance" associated with the
prediction. This quantity is the
variance of the difference between the "actual" value at that
location and our estimate of it - given our model of the parameter semi-variogram. This estimation variance may be related,
through its square root, to the traditional notion of "standard
error" of an estimation. In short, then, each estimate is accompanied by a second value measuring the potential
error or uncertainty at that location.
If we produce a grid of values, such as the one which
was used to construct Figure 2, then we also produce a grid of standard
errors. Thus a
second map may be constructed showing the uncertainty of the first estimated
map. Figure 3
shows the standard errors associated with the construction of the map of
potentiometric levels in Figure 2. We have used the term "standard
error" here to denote the square root of the estimation variance produced
during the Universal Kriging estimation process.
It can be seen from Figure 3
that the errors are lower in the area around the proposed repository, but rise
sharply away from this location. In
particular, they rise precipitously towards the
north-east-the very direction in which we
expect a freed particle to move. The
uncertainty associated with the travel path in this case rises almost
proportionally to its distance from the proposed repository. The solid line in Figure 2, then, can be seen to be merely a "best estimate" subject
to a great deal of potential error. We seek some way to relate the standard
error on the predicted potentiometric levels to the uncertainty on
the direction of travel of a radionuclide particle.
It has been
proposed that simulation may hold the key to the problem of quantifying
the uncertainty on the travel path. Two major approaches have been proposed. Several previous
papers (cf. Clifton & Neuman (1982), Dagan (1985), Devary et al. (l985)) suggests
that simulation of the fluctuations around the surface can be carried out using
the standard errors as shown in Figure 3. This approach would restore the
"stochastic" roughness of the surface to the smooth kriged grid and
allow a travel path to be calculated. If
many simulations are carried out, many possible travel
paths will result producing a "fan" of possible destinations for a
freed particle. This process requires a
simultaneous multivariate simulation of the errors at each point on the
estimation grid. A covariance matrix must be constructed between the errors at each node of the
grid and a realisation from this covariance matrix simulated. A direct implication of this approach is that
the errors at the grid points must be multivariate Normal. This approach may
also prove unwieldy in operation, since the simulation of a 12 by 12 grid of
points requires the simulation of a multivariate Normal from a 144 by 144 covariance matrix.
An alternative approach is to simulate the complete
surface at all grid points with the spatial structure of the original
parameters and the actual borehole values at those locations. This is known in geostatistics as 'conditional
simulation" and is well documented in the geostatistical literature
(cf. Journel and
Huijbregts (l978) and (1979)).
The underlying assumptions to geo-statistical conditional simulation are
of Normality for the original variables and of a particular model for the
semi-variogram of this data. In a case
such as the Pennsylvanian, where a strong trend exists, these assumptions apply
to the residuals from the trend surface.
Should the data be non-Normal, it is simple to a transformation on the
original data, construct the semi-variogram on this and then back-transform
once simulation has been carried out. A comparison of the two approaches to simulation is in
preparation (Harper &c Clark, 1989) and
will not be discussed here.
This paper discusses the latter "tradittional" conditional simulation
approach - that is, we simulate a surface which has the same spatial and statistical structure
as the original data and takes the actual data values at the actual sample
locations. The simulated surface will
have the same roughness as the original parameter data. If many surfaces are
simulated, an idea of the sensitivity of the travel path may be gained.
The Travel Path Fan
In addition to the "statistical" uncertainty
on the travel path of a freed radionuclide, we would like to discuss another
source of uncertainty which does not seem to have been
considered in the literature. The most
likely travel path has been illustrated in Figure 2 -
as a single solid line leading down the steepest slope on the potentiometric
surface. A simple intuitive analogy
would be that of water in a stream, which always follows the steepest path open
to it. However, this may be a simplistic
analogy. All we know
for certain is that flow will take place from a higher potentiometric
level to a lower one.
Then is no certainty that a particle will follow the
path of least resistance - it may flow to any lower level
which is available for it. Given
a fixed starting point, the path of a freed radionuclide may be any one of
several potential paths. Taking this
(admittedly simplistic) analogy, the likely travel path of a particle is in the
form of rather irregular fan spreading out from the original location. Figure 4 shows such a fan for the Kriged
potentiometric for the Pennsylvanian aquifer.
The contours represent the number of ways that a particular grid point
(on Figure 2) can be reached. That is, a value of '3' implies that this point is lower in
potentiometric level than three neighbouring points in potential travel paths.
The 'most likely' path has been superimposed
on this map. It can
easily be seen that the uncertainty on the most likely path is not,
symmetrical around it but is related to the dip on the estimated surface
itself.
This simplistic approach has been
based on the assumption that all paths downward are equally likely to be
taken by the particle. This is
fundamentally unrealistic and will be the subject of further research However,
for the purpose of a visual representation of the 'travel fan' this approach is
extremely useful since it shows all of the locations which could be reached by
a freed radionuclide particle.
Simulated
Travel Fans
In Figure 4, we saw that the most likely travel path of
a freed radio-nuclide particle actually covered an
area to the north-east of the proposed repository. The area which could be
rm6ed by a particle forms an irregular fan with its apex at the storage
site and following the general trend of the potentiometric levels as estimated
by the Universal Kriging method. This fan taken no account of the uncertainty inherent in the
estimation of the parameters at the grid points. This uncertainty has been
illustrated by mapping the standard errors for each estimated point
(Figure 3). For a more realistic idea of
the potential travel "fan" we must include
this uncertainty in the analysis.
Conditional simulation, in the geostatistical sense,
introduces the correct stochastic component back into the predicted surface and
shows a more realistic potentiometric surface - in terms of 'roughness'. This statement does not imply that a simulated surface
is a better estimator than a Kriged surface. This is simply one possible which honours the data values at the sampled
locations and has the same spatial and statistical structure as the original
data. This type of simulation has been
fully described elsewhere (cf Harper & Clark, 1988) and will be only
briefly here. Basically, such a surface is constructed
in three stages:
(i)
A Kriging
exercise is carried out on the original parameter data
to produce a grid of estimated points.
(ii)
Values are simulated at each sample location. The Kriging exercise is
repeated for the same grid.
(iii)
Values are simulated at each grid point.
The conditional simulation is produced
by subtracting (ii) from (iii) and adding the result to (i) for each grid point. Provided the simulated data followed the
semi-variogram and distribution as the original data, this will produce a
simulated surface which honours the data values and
has exactly the same characteristics as the original surface. In the case of a surface such as the
Pennsylvanian potentiometric data, the semi-variogram and distribution referred
to are than of the residuals from the trend surface.
In practice, the Kriging weights in stage (i) are stored
and used in stage (ii) to save computation time. In this way, the actual Kriging only has to be carried out once for any particular surface - no matter how many simulations are required. For each simulation, the travel path can be
calculated and accumulated.
The
Pennsylvanian Aquifer
The Kriged surface in Figure 2 was produced on a five mile grid.
Taking the same mesh size, 25 conditional simulations were
produced and the travel fans from each were stored. A study of individual simulations revealed
that some travel paths were extremely short with the particle becoming trapped in a local depression very soon after
leaving the repository. However, some
fans were much larger. Figure 5 shows
the accumulation of 25 such simulated travel fans.
Comparison with Figure 4 shows that the spread of the
simulated fan is much greater than that on the Kriged surface - as should be expected.
One interesting feature of the simulated fan is that there are possible
paths to the south and west of the repository.
On the kriged map these directions are
estimated to be "uphill". This
behaviour illustrates neatly one of the valuable contributions of the
simulation approach. Figure 2 shows the "best" estimation of the
potentiometric level at each one of the 5 mile grid
points. However, the true surface is
much more irregular than this smooth kriged surface. The simulation reminds us that there may be
directions of travel which would not even be
considered if only the estimated surface is considered.
Actual quantification of the risk associated with any
particular destination for a freed particle is a subject for more extensive
research than we have been able to cover in this study. However, we feel that this approach would
form a solid basis for developing the refinements required to determine
relative probabilities for all of the possible paths revealed by the present
study.
Practical
Aspects of Simulation
During the analysis which has been
described above, various factors were investigated as relevant to the
final results. In particular, the size
of the grid which was simulated appears to be a
definite influence on the size and shape of the final travel fan. The fan on the Kriged surface, of course,
remains the same - with minor refinements along
the outer envelope. However, the
simulated travel fans are heavily influenced by the size of mesh
which is studied. All of the
results above were produced on a five mile grid. Similar studies were also carried out for 2.5 mile grids and for a one mile grid. The fans tended to be similar in shape but
widely different as regards size. Figure
6 shows a comparison between the envelopes around 25 simulations on each of the
different size grids. At first sight
there seems to be a direct relationship between the refinement of the grid and
the area covered by the fan envelope. More detailed study reveals that the areas covered also
become more symmetrical around the repository site as the size of the grid
diminishes. This effect seems to be
rather puzzling until we recall the structure which is
being simulated.
The conditional simulation is producing a surface which is composed of three components: a very large
scale trend, dropping to the north-east; a Spherical semi-variogram stochastic
component, in which values are related up to 50 miles apart; and a completely
random
component measured by the nugget effect. At the five mile
scale, both trend and spherical components become reasonably strong so that the
travel paths are highly influenced by the trend and continuity on the simulated
surface. At the smaller scales, the nugget
effect becomes stronger than the continuity components resulting in more
"depressions" apparently trapping the freed particles. The one mile
simulation is almost symmetric around the site of the repository because the
simulated surface is almost entirely random - and therefore no preferred directions exist.
The major implications of this scale effect are that the
existing sample data is inadequate for this scale of investigation. There are very few borehole pairs within five
miles of one another. This means that
any geological variability on the one mile scale is
being lumped into the 'completely random' nugget effect. Quantification on this scale would require
much more detailed sampling to split this supposed nugget effect into its
proper spatial structures.
Multivariate
Conditional Simulation
We have seen that there are two
main components to the potential travel paths of a freed particle. The general
structure of the surface as revealed by the Kriging estimation gives the best
estimation of the potentiometric surface. This can be used
to produce a fan of possible travel paths available to a radionuclide. It may be possible to quantify the relative
probabilities of the various final destinations making up the envelope around
the travel path fan. The second
component is, paradoxically, the easier to quantify - being the stochastic variation which
has been smoothed out by the Kriging process, the uncertainty on the actual
estimated surface. This in measured by
the standard error on the estimation. More "realistic" surfaces may be simulated and the
deviation from the "best estimated" fan investigated.
We have discussed how that second component may be
quantified and simulated. The Kriging
process produces the smallest standard error, and hence the least amount of
stochastic variation to simulate. However, this last statement is only true if
we consider only the data available within the Pennsylvanian aquifer. It can be seen from Figure 1 that the north-east comer of the study area is particularly badly
sampled for values in the Pennsylvania aquifer. What we
cannot see from this Figure is that a companion aquifer, the Wolfcamp, which
lies above the Pennsylvanian is well sampled in this
area- In addition, there is a high correlation between the Potentiometric
levels in the Wolfcamp and those in the Pennsylvanian. An estimation method which
could take advantage of this extra information must (surely) be better
than a simple Universal Kriging on a single parameter.
An alternative to univariate Kriging is a multivariate
technique, generally known as Co-kriging.
With sample data such as that available in this
study traditional Co-Kriging cannot be applied.
An alternative approach, which has been dubbed Multivariate Universal
Co-Kriging (MUCK) has been expounded in detail in
Clark et al. (1988). It is sufficient to
say here that the technique is exactly analogous to Universal Kriging except
that the cross-correlation between samples in the Wolfcamp and values in the
Pennsylvanian may be included in the analysis.
The ultimate impact is to reduce the size of the estimation variances on
the Pennsylvanian estimates by a factor of up to 60%. The simulation process
which results includes a multivariate factor, since both Wolfcamp and
Pennsylvanian must be simulated, but in principle is similar to the univariate.
25 simulations were made on the MUCKed five mile grid on the Pennsylvanian aquifer. The travel path fans were
accumulated and we illustrated in Figure 7. Visual comparison of this envelope with that
in Figure 5 shows the reduction in uncertainty associated with the estimation
error on the potentiometric levels within the Pennsylvanian aquifer. This improvement more than compensates for
the minor increase in computer time associated with the MUCK estimation
method. In fact, if both aquifers were
to be simulated and studied in the investigation, the
computer time is considerably reduced over two completely separate studies.
Summary of
Discussion
In this paper we have
considered how geostatistical methods may aid in predicting the likely travel
path of a freed radionuclide from a high-level nuclear waste repository.
Traditional geostatistical estimation methods allow the prediction of the
'best' or most likely values for hydrological parameters at grid points over
the study area. From these the most likely travel path of a freed particle may
be calculated. In addition
a fan of possible paths may also be produced.
With the addition of geostatistical simulation,
fluctuations on the likely travel paths - as indicated by the smoothed Kriged
surface - may be studied. We have
illustrated this paper with only 25 simulations. Many more would be required for a definitive
study, but the principles are clearly illustrated. The possible destination of hazardous
material is seen to be a much larger area than would
be expected from a non-geostatistical study.
We have also shown that the uncertainty on the travel
path or fan can be reduced by introduction of other
parameters correlated to that one which is currently under study. Any number of associated parameters can be
included in such an investigation, provided that
sufficient data is available on each.
Whilst this case study has considered the release of radio-nuclide particles, it should be obvious that these
results are equally valid for any toxic or hazardous material which has been
stored close to a source of mobile groundwater.
References
Bair E.S., O'Donnell T.P. and Picking L.W., 1985,
"Hydrogeologic Investigation Based on Drill-Stem Test Data: Palo Duro
Basin Area, Texas and New Mexico", BMI/ONWI-566, prepared by Stone and
Webster Engineering Corporation for Office of Nuclear Waste Isolation, Battelle
Memorial Institute, Columbus, Ohio.
Clark I, Basinger K.Y. and Harper W.V., 1988, "MUCK
- A Novel Approach to Co-Kriging", Proc DOE/AECL '87 Conference on
Geostatistical, Sensitivity and Uncertainty Methods for Ground-Water Flow and
Radionuclide Transport Modelling, September, 1987.
Clifton P.M. and Neuman S.P., 1982, "Effects of
Kriging and Inverse Modelling on Conditional Simulation of the Avra Valley
Aquifer in Southern Arizona", Water Resources Research, Vol. 18,
No. 4, pp.l2l5-1234.
Dagan G., 1985, "Stochastic Modelling of
Groundwater Flow by Unconditional and Conditional Probabilities: The Inverse Problem",
Water Resources Research, Vol. 21, No. 1, pp.65-72.
Delhomme J.P., 1985, "Spatial Variability and
Uncertainty in Groundwater Flow Parameters: A Geostatistical Approach", Water
Resources Research, Vol. 15, No. 2, pp.269-280.
Devary J.L., Gupta S.I.K, Harper W.V., Cole C.R. and
Thompson F.L., 1985, "A Hybrid Approach to Uncertainty in Far-Field
Gw=dwater Flow", Proc. Symp. on Groundwater Flow and
Transport Modelling for Performance As t of Deep Geologic Disposal of
Radioactive Waste. A
Critical Evaluation of the State of the Art (HYDROCOIN), Albuquerque, NM, pp.233-256. .
Harper W.V., Basinger K.V. and Furr
J.M., 1986, "Geostatistical Analysis of Potentiometric Data in the
Pennsylvanian Aquifer of the Palo Duro Basin, Texas", BMI/ODWI-680, Office of
Nuclear Waste Isolation, Battelle Memorial Institute, Columbus, Ohio.
Harper W.V. and Clark I, 1988,
"Travel Path Uncertainty - A Case Study in Conditional Simulation",
Proc. of 3rd International Geostatistics Congress, Avignon, France.
Harper W.V. and Clark I., 1989, "Travel Path Uncertainty
by Geostatistical Simulation", in Proc. IGWMC
Conference on Groundwater Modelling, Indianopolis, February 1989.
Journel A.G. and Huijbregts Ch.J., 1979, Mining Geostatistics, Academic Press, New York.
