www.geologicacarpathica.sk
GEOLOGICA CARPATHICA, AUGUST 2009, 60, 4, 271—281 doi: 10.2478/v10096-009-0019-y
Stochastic hydrogeological modelling of fractured rocks:
a generic case study in the Mórágy Granite Formation
(South Hungary)
KÁLMÁN BENEDEK and GYULA DANKÓ
Golder Associates (Hungary) Ltd., Hűvösvölgyi út. 54, H-1021 Budapest, Hungary; kbenedek@golder.hu; gdanko@golder.hu
(Manuscript received June 3, 2008; accepted in revised form December 18, 2008)
Abstract: In connection with the Hungarian radioactive waste disposal program a detailed study of the mass properties of
the potential host rock (granite) has been carried out. Using the results of this study the various parameters (orientation,
length, intensity, transmissivity, etc.) describing a fracture set were estimated on the basis of statistical considerations.
These estimates served as basic input parameters for stochastic hydrogeological modelling of discrete fracture networks
(DFN), which is a strongly developing area of hydrology, providing geologically realistic geometry for site investigations.
The synthetic fracture systems generated were tested against some (but not all) field observations. The models built up on
the basis of the statistical descriptions showed the same equivalent hydraulic conductivity for the modelled region as the
field measurements. In addition, the models reproduce the observed hydraulic head-scattering along vertical boreholes. On
the basis of the stochastic simulations of the fracture system some input parameters for the performance assessment of the
planned repository were investigated. Calculation of flows into a planned disposal tunnel indicated that if the hydraulic
conductivity of the material in the tunnel is the only variable parameter then there are two thresholds: under 1
×10
—9
m/s
and above 1
×10
—5
m/s further change of the hydraulic conductivity does not dramatically affect the inflow.
Key words: Mecsek Mts, uncertainty, stochastic simulations, hydrogeological conditions, fractured rocks, granite.
Introduction
Fractured rock investigation often focuses mostly on discon-
tinuities. One reason, from a hydrogeological point of view,
is that fractures are usually more permeable than the sur-
rounding rock. Consequently flow and transport processes
are concentrated mainly in the fracture systems in the inves-
tigated rock mass. Fractures are also important in ore depos-
its, and petroleum, gas and geothermal reservoirs are
frequently located in fractured rocks (Van Golf-Racht 1982).
After the publication of the pioneering study by Snow
(1965), the hydrogeological approach to fractured rocks un-
derwent a gradual change. Researchers realized that fractures
behave as conducting features and often have a special role
in flow and transport processes.
The development of discrete fracture network (DFN) mod-
elling technologies was accelerated by underground research
laboratories (URL) for radioactive waste (Sweden, Canada,
United States, Japan, etc.). The aim was to characterize the
fractures and to develop tools able to reproduce measured
field data and to handle parameter uncertainties. These stud-
ies highlighted the need for proper integration of geological,
hydrogeological, geochemical and geophysical data (Neu-
man 2005). In addition, the hydrogeological studies in URLs
were able to provide some important input parameters for re-
pository performance assessments.
In most cases fractured rocks must be characterized in situ.
An advantage of in situ measurements is that the data in-
clude the effect of field stress conditions and of fracture con-
nectivity. The most widely used method is the single
borehole packer test and multi-borehole interference tests are
also used (Neuman 2005) in order to gain information about
the flow system on a larger scale. High-resolution flow-meters
(Rouhiainen & Pöllanen 1999) can provide detailed informa-
tion about the conductive fractures intersecting a borehole.
A key finding was that at some sites the fracture network
is not continuous (Uchida et al. 1998; Sawada et al. 2000).
This means that fractures or fracture clusters are not connect-
ed with each other, and thus that there is not hydraulic conti-
nuity between all points in a studied region. This feature is
known as compartmentalization.
Almost all the studies in URLs have concluded that all frac-
ture parameters (transmissivity, size, etc.) have a strong scale
dependence (La Pointe & Outters 2000; Andersson et al.
2002). For example, cm-scale fractures cannot be investigated
in regional scale mapping and, conversely, large tectonic lin-
eations cannot be observed in thin sections, although both of
them may belong to the same fracture system. These consider-
ations emphasize the need for mathematical techniques that
are able to integrate data measured at different scales.
Most approaches consider discontinuities in a fractured
rock mass as a two or three dimensional network of polygons
(Neuman 2005). These polygons are usually described by
stochastic parameters for location, size, orientation, trans-
missivity, etc., since these data cannot be determined for
each fracture deterministically (Dershowitz et al. 1998; M.
Tóth et al. 2004). The results from field measurements and
observations need some interpretation (conceptualization)
and statistical processing before becoming input parameters
for DFN modelling. The concept must be compatible with
the general geological and hydrogeological conditions of the
studied site.
272
BENEDEK and DANKÓ
The investigation of potential disposal sites for radioactive
waste in Hungary started in 1992 (Balla 2000). After a screen-
ing process the interest has been focused on Mississippian
granite (Mórágy Granite Formation) in the vicinity of the vil-
lage of Bátaapáti (Fig. 1). The geological exploration of the
area has been organized and financed by the Public Agency
for Radioactive Waste Management (PURAM). Several bore-
holes have been drilled (Balla et al. 2003) in the studied area
in order to obtain information about fracture orientation, frac-
ture density, rock stability etc. and to assess the geological/
hydrogeological suitability of the site (Balla et al. 2004). It is
important to note that this study is mainly generic and utilized
only data available from surface based investigations (Balla et
al. 2003) and did not consider data from access tunnel excava-
tions. The basis of this paper was originally presented by
Benedek & Mező (2005).
The structural evolution of the granite and the mechanisms
of different deformation phases are discussed by Maros et al.
(2004). The authors suggested that the pluton underwent a
complex tectonic evolution starting from the magmatic stage
and then subsequently the Variscan and the Alpine oroge-
nies. Two alternative theories were proposed to explain the
present structural pattern of the site: 1) folding related evolu-
tion; and 2) strike-slip faulting.
In the first section of this study we intend to reproduce the
geometry (orientation, length, density, etc.) of the fracture
system observed in the Mississippian granite formation and
additionally to estimate flow parameters of fractures (trans-
missivity). It is noteworthy that the up to date conceptual
model developed for the site will not be presented, since the
main aim of this paper is to demonstrate the capabilities of
the DFN approach. In the second step fracture systems are
built up stochastically (50 realizations) and we show differ-
ent, but not all aspects of the model’s applicability.
Site geology
The site discussed in this paper is located in the Mecsek Mts
in southern Hungary (Fig. 1). On the basis of geological inter-
pretations the system can be divided into two main segments:
1) overlying sediments and 2) the granite system (Balla et al.
2003). It is important to note that the granite can also be divid-
ed into subsystems: 1) weathered granite and 2) fresh granite.
This paper focuses on the latter one.
Buda et al. (2000) distinguished four different types of crys-
talline rocks at the investigated site: 1) microcline granitoids
containing megacrysts, 2) amphibole-rich enclaves, 3) micro-
granites, 4) pegmatites. They suggested, based on whole-rock
major and trace element geochemistry, that partial melts were
formed in a continental collision zone during the Variscan
orogeny. In addition, basaltic melts could have been generated
as a consequence of partial fusion of a mantle wedge above
the subduction zone. The granite body underwent numerous
secondary processes: K-metasomatism, mylonitization, cata-
clasis, hydrothermal alteration and mineralization. In addition,
the granite was affected by multiphase metamorphic events
(Koroknai 2003) which resulted in the formation of slate, my-
lonite, etc. The age of the granite formation is the subject of
debate: Balogh et al. (1983) suggest a Carboniferous forma-
tion age, however Chernysev et al. (2002) argues for a meta-
morphic event at that time.
The granite was penetrated by swarms of Cretaceous ba-
saltic (bostonite, trachyandesite) dykes (Harangi 2003). At
Fig. 1. Generalized map of the study area.
273
STOCHASTIC HYDROGEOLOGICAL MODELLING OF FRACTURED ROCKS (SOUTH HUNGARY)
the interface between the granite and dykes a cooling margin
can frequently be observed and this is usually associated
with fractures. The dykes are strongly weathered and hydro-
thermally altered.
Hydrogeological conditions
At the site of the planned repository very detailed investi-
gations have been carried out using various methods. The
hydrogeological conditions of the entire investigated area is
presented in the basic paper of Balla et al. (2004).
Among the features of the fresh granite, one of the most im-
portant is its strongly fractured and argillized character in local
zones of tectonic displacements. It is also heterogeneous (Ma-
ros et al. 2003). Extensive single-hole packer testing, with an
average test interval length of 10 m, found that packer interval
transmissivities vary over six orders of magnitude without any
spatial trends (Balla et al. 2004). The calculated transmissivity
of the fresh granite is in the range between 10
—12
and 10
—6
m
2
/s.
The transmissivity of the main water conductive features var-
ies between 8
×10
—6
and 2
×10
—5
m
2
/s (Balla et al. 2004).
A typical hydraulic head profile along a borehole is dis-
played in Fig. 2. The variability of hydraulic head results
from two different phenomena: 1) small scale (1 to 5 m) hy-
draulic head scattering and 2) large (5 to 20 m) hydraulic
head jumps. The first is caused by the natural variability
(fracture orientation, size, transmissivity and connectivity)
and clustering of fracture systems. Balla et al. (2004) suggest
that the large hydraulic head jumps are associated with the
strongly altered, argillized fault core zones observed in the
granite. It is frequently concluded that the boundary between
compartments with different hydraulic heads must have a
very low permeability (Dershowitz et al. 1998). In the case
of the studied area Tóth et al. (2003) estimated the hydraulic
conductivity of these zones to be 10
—12
to 10
—10
m/s by using
inverse modelling technique. Cross-hole interference tests
clearly demonstrated that borehole sections displaying only
small scale hydraulic head-scattering are within a single
compartment (Benedek et al. 2003a). At site scale the strong-
ly compartmented character causes high complexity of the
flow pattern. In this paper the effects of compartment bound-
aries and the general conceptualization of the site are not
considered further. This paper studies conditions within a
single compartment.
The chemistry of groundwater at the site has been studied
in detail by Horváth et al. (2003). They suggested that in the
northern area the vertical flow velocity component is larger
than in the South. This is presumed to reflect differences in
the compartmentalization of the granite. In this paper only
boreholes located in the southern section were considered. In
this part of the site 6 boreholes (Üh-3, Üh-4, Üh-5, Üh-22,
Üh-26, Üh-28) have been drilled (Fig. 1).
Modelling tool applied
The modelling tool (FracMan software package) applied
was developed by Golder Associates Inc. (Seattle) to model
the geometry of discrete features, including faults, fractures,
paleochannels, karstic elements and stratigraphic contacts
(Dershowitz et al. 1998). This software is frequently used in
nuclear projects and in the oil industry. The oil industry ap-
plies the software in up-scaling information from fractured
rock reservoirs to equivalent continuum models, to predict
well production and to provide information about the neces-
sary well spacing.
FracMan simulates fractures as 3-D networks of triangular
finite elements. The software represents individual fractures
as 2-D elements defined stochastically based on probability
distributions for orientation, length, transmissivity, aperture,
storativity, etc. In addition, deterministic features can also be
involved in the models. The solver for FracMan (MAFIC)
uses a Galerkin finite element solution scheme to approxi-
mate the solution for the diffusivity equation (Bear 1972).
FracMan is also able to model coexisting matrix (generic
matrix block or fully discretized matrix elements) and frac-
ture (DFN elements) systems. Solute transport modelling un-
der steady-state or transient flow conditions, and stochastic
emulation of convective dispersion, matrix diffusion, radio-
active decay and mineral specific retardation are also sup-
ported. Boundary conditions are assigned in the traditional
way (constant or transient hydraulic head and flux).
FracMan provides an integrated environment for the entire
process of discrete feature data analysis and modelling. In
Fig. 2. Head profile along a borehole (Üh-28) at the site. Small scale
head-scattering and larger jumps in head are indicated with arrows.
274
BENEDEK and DANKÓ
the first step the user can analyse and transform raw field
data (fracture orientation, packer interval transmissivity,
fracture density, etc.) into the format required for discrete
fracture modelling. In the next step, stochastic simulations of
fracture patterns can be made. If the generated fracture sys-
tems meet the selection criteria they can be converted into a fi-
nite element mesh with the appropriate hydraulic parameters.
The software is also able to model fractures generated un-
der specific geological and geodynamic conditions. For ex-
ample, the software has the capability to reproduce fractures
generated due to folding, faulting, hydrofracturing or result-
ing from a specific in situ stress field.
Model size, boundary conditions
The size and boundary conditions for the model were
highly simplified, recognizing that for a specific location
compartment boundaries may need to be considered. The
model was a generic cube with 200 m sides. The boundary
conditions were defined to result in a one directional flow of
1% gradient, from S to N. The other faces of the cube were
defined as no-flow boundaries.
The model contained stochastic and deterministic ele-
ments. The stochastic elements represent the fractures in the
granite and the deterministic elements the planned waste dis-
posal tunnel.
Stochastic elements
Fracture orientation
The total set of fracture orientation data was measured in
individual boreholes by BHTV (acoustic borehole teleview-
er; Zilahi-Sebess et al. 2003). For this study only the hydrau-
lically active fractures have been considered, as other ones
have only a limited role in the flow system. The selection
was made using the results of high-sensitivity flow measure-
ments (Heat-Pulse Flowmetry (HPF); Szongoth & Galsa
2003). These were carried out at 1 m scale, and provide data
on inflow intensity and location. The combination of the two
data sets (BHTV and HPF) allows the estimation of orienta-
tions and of the one-dimensional (along boreholes) fracture
intensities of the hydraulically active fractures.
Visually two (or three) sets of fractures can be distinguished:
1) NE-SW strike, SE dip; 2) E-W strike, N dip; 3) N-S, NW-SE
strike, W, SW dip (Fig. 3). For this study the probability dis-
tributions of the fracture orientation data were approximated
using the univariate Fisher distribution (Mardia 1972). The
probability distribution is defined as follows:
where
κ is the distribution parameter and ϕ’, Θ’ are the
variations about the mean direction. The goodness of fit was
demonstrated by Kolmogorov-Smirnov statistics. To reduce
the overestimation of horizontal and sub-horizontal fractures
due to the vertical orientation of the boreholes the Terzaghi
(1965) correction was applied.
For the study described in this paper only two sets of frac-
tures (NE-SW and N-S, NW-SE strike) were defined in order
to simplify the system and to reduce computation time. The
results are presented in Table 1.
Fracture size distribution
The 3-D size distribution of the fractures in the rock mass
can be estimated by different methods: 1) size-number scal-
ing plots at regional, outcrop, etc. scales; 2) censoring statis-
tics; 3) trace length matching; 4) pressure transient analysis.
In this study the size-number scaling plot approach was used.
At the site and in the vicinity of the site there are three in-
dependent trace length data sets: 1) artificial outcrops (Balla
et al. 2003; Gyalog et al. 2003); 2) vertical seismic profiling
(VSP) in boreholes (Prónay 2003); 3) seismic reflection sur-
veys (Prónay et al. 2003).
Fig. 3. Stereographic projection (lower hemisphere) of hydraulical-
ly active fractures in the boreholes investigated.
Fisher parameters
Fracture system
Strength (%)
Main dip direction (
°)
Main dip angle (
°)
Fisher K
Kolmogorov-Smirnov
statistics (%)
Set1 78.6
341.6 82.5
5.71 95.6
Set2 21.4
126.1 59.5
27.32 93.6
Table 1: The statistical parameters (Fisher distribution) for the orientation of the two hydraulically active fracture sets.
(
)
π
π
ϕ
κ
ϕ
κ
ϕ
κ
2
'
0
,
)
1
(
2
'
sin
'
,'
'
cos
≤
Θ
≤
−
=
Θ
e
e
f
275
STOCHASTIC HYDROGEOLOGICAL MODELLING OF FRACTURED ROCKS (SOUTH HUNGARY)
These data sets are representative of different scales: out-
crops provided information about the trace length up to
~ 6 m, VSP in the range of 2 to 40 m, and seismic surveys in
the range of 50 to 380 m. Trace length distributions at any
scale are truncated in two ways: 1) the lower one is selected
by the operator; 2) the upper one results from the size of the
sampling window.
The scaling properties of trace length have been estimated
using the complementary cumulative density function
(CCDF; Fig. 4) for each scale. Due to the truncations only
the linear segments of these functions can be applied in the
estimation of fracture size distributions (Andersson et al.
2002). At any scale the best-fit to the measured data was a
power law function (Fig. 4). The data for different scales are
very consistent, which suggests that the fractures at all scales
belong to the same system.
Direct simulation of fractures in 3-D based on 2-D data
sets is not reliable. For example, the probability that a large
fracture intersects an outcrop is higher than is the case for
small fractures, and consequently trace length is biased. In
addition, the orientation of fractures and outcrops can also
influence the trace length distribution. La Pointe & Outters
(2000) argued for a relationship between the size distribution
of trace length (2-D) and that of fractures (3-D) intersecting a
trace plane (outcrop). The following equations can be ap-
plied to link the parameters of fracture size distribution and
trace length distribution:
D
3-D, fractures
= D
2-D, traces
+ 1,
x
3-D, fractures
= x
2-D, traces
where D
3-D, fractures
and x
3-D, fractures
are the parameters of the
power law function in 3-D and D
2-D, traces
and x
2-D, traces
are the
parameters of a power law function in 2-D. D and x represent
the exponent and minimum value of power law distribution,
respectively. In this study D was calculated as an average of D
values at different scales (D = 3.26) but x was determined as
the smallest value among the different scales (x = 0.145).
The size distribution of fractures with different orientation
was also compared. This complex data set (orientation, size) is
only available for artificial outcrops (Balla et al. 2003; Gyalog
et al. 2003). CCDF calculations suggested that there is no fun-
damental difference between the trace lengths of differently
oriented fractures (Benedek & Mező 2005).
Fracture intensity
The intensity of fractures in 3-D can be described as a vol-
umetric value, as fracture area per unit volume of rock (P
32
,
m
2
/m
3
). However this parameter cannot be estimated direct-
ly from field measurements and must be derived indirectly.
The number of fractures intersecting individual boreholes
(linear intensity, fractures per meter, P
10
) is only a function
of the orientation of boreholes, fractures and the size of frac-
tures (La Pointe et al. 1993). Dershowitz et al. (1998) sug-
gested that volumetric and linear intensity can be linked
through the following equation:
P
32
= P
10, obs
× (P
32, sim
/P
10, sim
) = C
31
× P
10, obs
C
31
= P
32, sim
/P
10, sim
where P
32, sim
is the volumetric intensity in the simulation,
P
10, obs
is the linear intensity measured in the boreholes (in
this study BHTV data are used; Zilahi-Sebess et al. 2003)
and P
10, sim
is the linear intensity in the simulation. C
31
is a
proportional constant, which can be determined by using in-
verse modelling. Since the size distribution of the fractures
and the orientation of the boreholes and the fractures are
known, all the data required to calculate P
32
are available.
P
32, sim
was changed gradually until an acceptable fit between
Fig. 4. The complementary cumulative density function (CCDF) of trace length for different scales.
276
BENEDEK and DANKÓ
P
10, obs
and P
10, sim
was found. In the stochastic simulations
with 50 realizations the average value of C
31
in the individu-
al realizations stabilized at around 1.35 for both fracture ori-
entation sets (Benedek & Mező 2005).
In general only a portion of the entire fracture system is
active hydraulically, since the majority are very small (back-
ground fracturing) or not connected to other fractures. For
this reason, and to save computational time, the size distribu-
tion was truncated. The upper limit was determined on the
basis of the cross-hole interference tests, which indicate that
the largest distance between packer intervals in hydraulic
connection is about ~ 400 m (Benedek et al. 2003b). The
lower limit was estimated using inverse modelling: the lower
limit for fracture size was changed until the modelled linear
intensity of hydraulically active fractures was consistent
with field measurements in boreholes (HPF measurements in
Szongoth & Galsa 2003). Based on this the lower limit value
was set to 10 m. The comparison of linear intensity indicated
by HPF measurements and the simulated values in 50 indi-
vidual realizations is shown in Fig. 5.
Spatial model
The spatial model is basically the definition of the rela-
tionship of the individual fracture centres in space, for which
many models are available (Dershowitz et al. 1998). Bene-
dek et al. (2003b) carried out a series of fractal dimension
calculations for each borehole in order to evaluate the spatial
relationships. On the basis of the location and the frequency
of fracture-borehole intersections, box-fractal and mass-frac-
tal dimensions were calculated. These calculations showed
fractal dimensions (Mandelbrot 1985) in the range of 0.88 to
0.99, very close to one. A fractal dimension of 1 in a 1-D
data source (borehole) represents a random fracture pattern
in space. Therefore a Baecher model (Baecher et al. 1977), a
randomly, but uniformly distributed fracture centre model
was applied.
Transmissivity distribution of individual fractures
For a DFN flow simulation it is necessary to estimate the
transmissivity distribution of individual fractures. For this
study an approach adapted from Osnes et al. (1988) was
used. This method assumes that the transmissivity of a pack-
er test interval can be regarded as the sum of the transmissiv-
ities of the conductive fractures intersecting the test interval:
∑
j
n
1
j=
j
i,
i
T
=
T
where T
i
is the transmissivity of the packer interval, n
j
is
the number of conductive fractures intersecting the test inter-
val and T
i,j
are the transmissivities of the conductive frac-
tures. The approach used was the following: 1) conductive
fractures with an intensity approximated by a Poisson distri-
bution are generated stochastically for the test zone, 2) to the
conductive fractures transmissivities are assigned stochasti-
cally by Monte-Carlo simulation based on an initial estimate
of the fracture transmissivity distribution (in this study a log-
normal distribution), and thus the interval transmissivity can
be calculated, 3) stochastic simulation terminates when a
good fit is found between observed packer interval transmis-
sivities and simulated transmissivities. The goodness-of-fit
is tested using Kolmogorov-Smirnov statistics.
Test interval transmissivities are available from three inde-
pendent sources: short-term single-hole tests (Molnár et al.
2000, 2003a,b), 2) long-term single-hole tests (Molnár et al.
2000, 2003a,b) and 3) cross-hole interference tests (Bradley
et al. 2000; Ács et al. 2003a,b,c). The relatively small num-
ber of interference tests precludes a reliable estimate of the
transmissivities and therefore only the single-hole tests were
used. In the simulations the following assumptions were ap-
plied: 1) the P
10
value (linear intensity) of conductive frac-
tures is close to 0.08 based on HPF measurements (Szongoth
& Galsa 2003), 2) the detection limit of HPF measurements
Fig. 5. Comparison of field
P
10
measurements and simu-
lated values from individual
realizations with fracture
length distribution truncated
below 10 m and above
400 m. Average – stabili-
zation of average values with
realizations; realization –
calculated value of individu-
al realizations; measured –
field measurements.
277
STOCHASTIC HYDROGEOLOGICAL MODELLING OF FRACTURED ROCKS (SOUTH HUNGARY)
was ~ 1
×10
—8
m
2
/s and therefore test intervals with transmis-
sivity below this threshold were considered non-conductive.
The results of simulations for short-term and long-term tests
are presented in Table 2. In the further modelling we used
lognormal transmissivity parameters derived from long-term
measurements.
Deterministic elements
For this study some deterministic continuum elements
were generated to represent a future waste disposal tunnel.
The tunnel was placed in the centre of the modelled block
so that groundwater flows in the direction perpendicular to
it. The longitudinal axis of the tunnel was 100 m long and its
cross-section was assumed to be a square with 10 m sides.
The hydraulic behaviour of the tunnel content was described
by equivalent porous continuum elements (hydraulic con-
ductivity varied from 1
×10
—12
m/s to 1
×10
—3
m/s) with a cell
size of 2.5 m
×2.5 m×2.5 m.
Results
On the basis of the geometric and flow parameter distribu-
tions 50 stochastic realizations were generated. The fracture
system for one realization is presented in Fig. 6. Different
aspects of the 50 realizations are discussed by Benedek &
Mező (2005), in this paper only some selected results are
presented.
Equivalent hydraulic conductivity of the modelled region
The equivalent hydraulic conductivity of the modelled re-
gion was calculated on the basis of the flow through the con-
stant hydraulic head boundaries. Since the hydraulic gradient
was set to be 1%, the application of Darcy’s law allows the
calculation of the equivalent hydraulic conductivity.
Balla et al. (2004) suggested that the hydraulic conductivi-
ty of the fractured granite varies between 1
×10
—12
m/s and
1
×10
—6
m/s at the scale of 10 m packer intervals and that it
has a lognormal distribution with a mean of 3
×10
—9
m/s. The
Table 2: Statistical parameters for individual fracture transmissivities derived using the method from Osnes et al. (1988).
Short-term tests
Long-term tests
Mean (m
2
/s), µ
logT
3.2×10
–7
8.0×10
–7
Deviation (m
2
/s), σ
logT
4.2×10
–6
3.4×10
–6
Fracture intensity (P
10
)
0.0881 0.0891
The ratio of fractures with a transmissivity less than 1×10
–8
m
2
/s
52.6 % in the model; 47.4 % in field
measurements
16.1 % in the model; 9.68 % in field
measurements
Kolmogorov-Smirnov test
98.4 95.9
Fig. 6. Calculated head distribution in a single realization.
278
BENEDEK and DANKÓ
median transmissivity value for the boreholes considered in
this study (southern section of the site, Fig. 1) is 1.4
×10
—9
m/s.
The calculated results for the 50 stochastic realizations indicat-
ed a median equivalent hydraulic conductivity of 6
×10
—9
m/s,
and thus it can be concluded that the model matches the field
measurements well.
Flow patterns calculated
Field studies (Ács et al. 2003d) indicated that one of the most
characteristic hydrogeological features of the site is the fre-
quently observed smaller-scale hydraulic head-scattering and
the larger hydraulic head jumps along individual boreholes
(Fig. 2). In this paper the larger hydraulic head jumps and their
conceptual origin (in compartmentalization) are not considered.
Figure 6 shows results for a single realization. The calcu-
lated hydraulic head does not change gradually from the low
hydraulic head boundary to the high hydraulic head one, but
it is strongly affected by system heterogeneity. In order to in-
vestigate the hydraulic head heterogeneity nine virtual verti-
cal boreholes were inserted and hydraulic head profiles were
calculated along them (Fig. 7). The hydraulic head profile
for borehole 6 shows a very high head difference (from
~100 m to ~98 m), which covers almost the entire range be-
tween the boundary conditions. On the other hand, borehole
3 has a very smooth hydraulic head profile. The interpreta-
tion of different hydraulic head profile patterns is that a clus-
Fig. 7. Head profiles along some vertical boreholes. The numbers
refer to the boreholes inserted into the modelled region.
Fig. 8. Average, minimum and maximum values of inflow through
the walls of the disposal tunnel as a function of tunnel content con-
ductivity.
ter of fractures may be connected to the high hydraulic head
boundary or to the low hydraulic head boundary, and that
different fractures or fracture clusters are not necessarily
connected directly, but only through the rock matrix. This
observation emphasizes the importance of connectivity and
geometry in fractured rocks. This conclusion also demon-
strates that a fractured rock mass cannot be characterized by
a constant, uniform hydraulic gradient, as the actual gradi-
ents are highly scale and location dependent.
Flow into the tunnel
One of the most important parameters for performance as-
sessment obtainable from hydrogeological models is the
flow through the disposal tunnels, as groundwater can trans-
port radionuclides to the geosphere and biosphere. The flow
is influenced by three engineering elements: 1) the material
in the tunnel, 2) the excavation disturbed zone, 3) grouting
of transmissive fractures. In the study described here only the
effect of the material in the tunnel was investigated. The other
two elements were studied by Benedek & Mező (2005).
The calculated average, minimum and maximum inflow
values through the faces of the drift are displayed as a func-
tion of tunnel content conductivity in Fig. 8. There is a
threshold under which the reduction of the hydraulic conduc-
tivity does not significantly reduce the flow. This threshold
is ~ 1
×10
—9
m/s. There is also an upper threshold, above
which increasing the hydraulic conductivity does not in-
crease the flow further. This threshold is ~ 1
×10
—5
m/s.
Flux inside a tunnel
The flux distribution within a disposal tunnel may be sig-
nificant for performance assessment, since in a fractured
rock mass concentrated inflow and outflow could occur
close to fracture-tunnel intersections and result in heteroge-
neous conditions. The drift was divided into 40 sections
along the longitudinal axis and the flux through each drift
section was calculated.
The calculated flux distribution along the longitudinal axis
is shown in Fig. 9A,B for one of the 50 realizations. This fig-
279
STOCHASTIC HYDROGEOLOGICAL MODELLING OF FRACTURED ROCKS (SOUTH HUNGARY)
Fig. 9. The x-direction flux in the 40 sections of the drift based on one of the 50 realizations in the cases of low (A – 1
×10
—12
m/s) and
high (B – 1
×10
—3
m/s) backfill conductivity. The negative or positive signs of the values calculated refer to the direction of flux.
280
BENEDEK and DANKÓ
ure clearly demonstrates that the flux distribution is not uni-
form if the tunnel content conductivity is low (1
×10
—12
m/s):
flux intensity is much higher close to fracture-drift intersec-
tions than far away from intersections. In the case of low hy-
draulic conductivity (Fig. 9A) some abrupt changes between
neighbouring sections can be observed, however, in the case
of high hydraulic conductivity (Fig. 9B) the distribution along
the axis is more smoothed. In that case the tunnel content con-
nects individual fracture-tunnel intersections and significant
longitudinal flux is generated. It may be noted that the small
flux values in Fig. 9A,B suggest that the most important trans-
port process inside the tunnel is expected to be diffusion.
Conclusions
To assess a repository for radioactive waste it is crucial to
have a good estimate of the flow of groundwater through the
waste, since groundwater is able to transport radionuclides
into the geosphere and biosphere.
In the study described here, based on the available geolog-
ical and hydrogeological data sets, the most important pa-
rameters describing the hydraulically active fracture systems
in the host rock were estimated (orientation, intensity,
length, transmissivity, etc.) on a statistical basis. The fracture
system was then approximated using a DFN model. In the
model the waste disposal tunnel was represented using
equivalent porous continuum elements.
In this paper only a few aspects of the DFN model and a
selection of results have been displayed in order to demon-
strate the capabilities of this technology.
The most important conclusions of the work can be sum-
marized as follows:
The hydraulically active fractures have a NE-SW strike
and the orientation can be described very well with two Fisher
distributions. The directions coincide with the main tectonic
orientations observed in the granite.
The fracture length distribution can be characterized
with a truncated power-law function based on the trace
length measurements at different scales.
The transmissivity distribution of individual fractures
can be defined by a lognormal distribution based on the re-
sults of short-term and long-term single-borehole tests.
The equivalent hydraulic conductivity based on field
measurements was well reproduced by the model.
The model displayed the observed small-scale stochastic
hydraulic head-scattering along boreholes.
The inflow is determined mainly by the content of the
emplacement tunnel. The larger its hydraulic conductivity,
the larger the flow into the tunnel.
The flux distribution inside the tunnel is heterogeneous
to a degree depending on the hydraulic conductivity of the
tunnel content. The smaller the hydraulic conductivity the
more pronounced the heterogeneity.
The flux inside the tunnel is very low, which means that
diffusion is likely to be the most important transport process.
Acknowledgments: The authors would like to thank Zoltán
Bőthi, Péter Molnár, Gyula Mező (Golder Associates Hungary
Ltd.) for their constructive comments in the preparation period
of the manuscript. We are also grateful to Martin Goldsworthy
(Golder Associates GmbH, Celle) for strong improvement of
the manuscript’s language level. This paper could be pub-
lished with the permission of the Public Agency for Radioac-
tive Waste Management (PURAM). We also wish to thank the
three reviewers and journal editors, whose comments and sug-
gestions substantially improved this manuscript.
References
Andersson J., Berglund J., Follin S., Hakami E., Halvarson J., Her-
manson J., Laaksoharju M., Rhén I. & Wahlgren C. 2002:
Testing the methodology for site descriptive modelling. Swed-
ish Nuclear Fuel and Waste Management Co., Technical Re-
port, TR-02—19.
Ács V., Molnár P. & Enachescu C. 2003a: Documentation of the
southern interference test. Manuscript, Bátatom Ltd., Budapest,
BA-03—219 (in Hungarian).
Ács V., Molnár P. & Enachescu C. 2003b: Documentation of the
northern interference test. Manuscript, Bátatom Ltd., Budapest,
BA-03—220 (in Hungarian).
Ács V., Molnár P. & Enachescu C. 2003c: Documentation of the
middle interference test. Manuscript, Bátatom Ltd., Budapest,
BA-03—221 (in Hungarian).
Ács V., Benedek K., Mező Gy. & Molnár P. 2003d: Integrated hydro-
geological interpretation. Hydraulic potentials at the site. Manu-
script, Bátatom Ltd., Budapest, BA-03—216 (in Hungarian).
Baecher G.B., Lanney N.A. & Einstein H.H. 1977: Statistical de-
scription of rock properties and sampling. Proceedings of the
18th U.S. Symposium on Rock Mechanics. Amer. Inst. Mining
Engineers 5C1—8.
Balla Z. 2000: Exploration and characteristics of the Üveghuta site.
Ann. Report Geol. Inst. Hung., 59—77.
Balla Z., Albert G., Chikán G., Dudko A., Fodor L., Forián-Szabó
M., Földvári M., Gyalog L., Havas G., Horváth I., Jámbor Á.,
Kaiser M., Koloszár L., Koroknai B., Kovács-Pálffy P., Maros
Gy., Marsi I., Palotás K., Peregi Zs., Rálisch L.-né, Rotárné
Szalkai Á., Szőcs T., Tóth Gy., Turczi G., Prónay Zs., Vértesy
L., Zilahy-Sebess L., Galsa A., Szongoth G., Mező Gy., Mol-
nár P., Székely F., Hámos G., Szűcs I., Turger Z., Balogh J.,
Jakab G. & Szalai Z. 2003: Final report of surface based inves-
tigation, Bátaapáti (Üveghuta), 2002—2003. Manuscript, Báta-
tom Ltd., Budapest, BA-03—156 (in Hungarian).
Balla Z., Horváth I., Tóth Gy., Benedek K., Mező Gy. & Molnár P.
2004: Hydrogeological pattern of the Bátaapáti (Üveghuta)
site. Ann. Report Geol. Inst. Hung., 449—472.
Balogh K., Árva-Sós E. & Buda Gy. 1983: Chronology of granitoid
and metamorphic rocks of Transdanubia (Hungary). An. Inst.
Geol. Geofiz. 61, 359—364.
Bear J. 1972: Dynamics of fluids in porous media. Amer. Elsevier
Publ. Co., New York, 1—764.
Benedek K. & Mező G. 2005: Hydrogeologic modelling. Flow cal-
culations through deterministic surfaces. Manuscript, Bátatom
Ltd., Budapest, RHK-K—055/05 (in Hungarian).
Benedek K., Ács V., Andrássy M. & Molnár P. 2003a: Integrated
hydrogeological interpretation. Hydraulic connections at the
site. Manuscript, Bátatom Ltd., Budapest, BA-03—218 (in Hun-
garian).
Benedek K., Outters N. & Hermanson J. 2003b: Hydrogeological
modelling for performance assesment. FracMan model. Manu-
script, Bátatom Ltd., Budapest, BA-03—178 (in Hungarian).
Bradley J.G., Enachescu C., Macdonald B. & Molnár P. 2000: Hy-
281
STOCHASTIC HYDROGEOLOGICAL MODELLING OF FRACTURED ROCKS (SOUTH HUNGARY)
drogeological interference testing in the Carboniferous granite,
at Üveghuta, South-west Hungary. Ann. Report. Geol. Inst.
Hung., 427—437.
Buda Gy., Puskás Z., Gál-Sólymos K., Klötzli U. & Cousens B.L.
2000: Mineralogical, petrological and geochemical characteris-
tics of crystalline rocks of the Üveghuta boreholes (Mórágy
Hills, South Hungary). Ann. Report Geol. Inst. Hung., 231—242.
Chernysev I., Volkov V., Mohov V., Lapina A., Jakusev M., Dubiny-
in A., Kogan A., Lebegyev Sz., Arakeljanc V., Nyikisina M. &
Satagin K. 2002: K-Ar and Rb-Sr geochronology of the Mórágy
granite. Manuscript, Bátatom Ltd., Budapest, BA-02—56.
Dershowitz W.S., Lee G., Geier J.E., Foxford T., La Pointe P. &
Thomas A. 1998: FracMan. Interactive discrete feature data
analysis, geometric modelling, and exploration simulation.
User documentation. Version 2.6. Golder Associates Inc., Se-
attle, Washington, 1—189.
Gyalog L., Jámbor Á., Kókai A., Maros Gy., Peregi Zs., Konrád Gy.,
Máthé Z. & Szebényi G. 2003: Geological description of A1 and
A2 outcrops. Manuscript, Bátatom Ltd., Budapest, BA-03—78.
Harangi Sz. 2003: The petrography of rock samples from the Me-
csek Mts. and their comparison with E-Mecsek Lower Creta-
ceous volcanic rocks. Manuscript, Hung. Geol. Inst., Budapest.
Horváth I., Marsó K., Muráti J., Nagy P., Rotárné Szalkai Á., Szőcs
T. & Tóth Gy. 2003: Integrated hydrogeological interpretation.
Manuscript, Bátatom Ltd., Budapest, BA-03—123.
Koroknai B. 2003: Microtectonic investigation of oriented rock
samples. Manuscript, Bátatom Ltd., Budapest, BA-03—84.
La Pointe P. & Outters N. 2000: Evaluation of the conservativeness
of the methodology for estimating earthquake-induced move-
ments of fractures intersecting canisters. Swedish Nuclear Fuel
and Waste Management Co., Technical Report TR-00—08.
La Pointe P.R., Wallmann P.C. & Dershowitz W.S. 1993: Stochas-
tic estimation of fracture size from simulated sampling. Int. J.
Rock Mechanics, Mining Sci., Geomechanics Abstr., Vol. 30,
1611—1617.
M. Tóth T., Szűcs É., Schubert F. & Hollós Cs. 2004: Conceptual
fracture network model of the crystalline basement of the Sze-
ghalom Dome (Pannonian Basin, SE Hungary). Acta Geol.
Hung. 47/1, 19—34.
Mandelbrot B.B. 1985: Self-affine fractals and fractal dimension.
Physica Scripta 32, 257—260.
Mardia K.V. 1972: Statistics of directional data. Academic Press,
New York, 1—357.
Maros G., Balla Z., Dudkó A., Fodor L., Forián-Szabó M., Koroknai
B., Lantos M. & Palotás K. 2003: Final report: Tectonics.
Kézirat, MÁFI, Budapest, Tekt. 1046; Manuscript, Bátatom
Ltd., Budapest, BA-03—118.
Maros G., Koroknai B., Palotás K., Fodor L., Dudko A., Forián-
Szabó M., Zilahi-Sebess L. & Bán-György E. 2004: Tectonic
and structural evolution of the north-eastern Mórágy Block.
Ann. Report Geol. Inst. Hung., 2003, 370—386.
Molnár P., Bradley J.G., Enachescu C. & Wozniewicz J. 2000: Sin-
gle-borehole hydrogeological testing in the Carboniferous gran-
ites, at Üveghuta, in South-west Hungary. Ann. Report Geol.
Inst. Hung., 407—417.
Molnár P., Ács V., Andrássy M., Rőczei N., Szűcs N. & Enachescu
C. 2003a: Final report on the hydraulic testing of well Üh—26.
Manuscript, Bátatom Ltd., Budapest, BA-03—201.
Molnár P., Ács V., Andrássy M., Rőczei N., Szűcs N. & Enachescu
C. 2003b: Final report on the hydraulic testing of well Üh—28.
Manuscript, Bátatom Ltd., Budapest, BA-03—83.
Neuman S.P. 2005: Trends, prospects and challenges in quantifying
flow and transport through fractured rocks. Hydrogeol. J. 13,
124—147.
Osnes J.D., Winberg A. & Andersson J. 1988: Analysis of well test
data – application of probabilistic models to infer hydraulic
properties of fractures. Topical Report RSI-0338, RE/SPEC
Inc., Rapid City, South Dakota.
Prónay Zs. 2003: Report on radar measurements in boreholes.
Manuscript, Bátatom Ltd., Budapest, BA-03—37A.
Prónay Zs., Neducza B. & Törös E. 2003: P- and S-wave seizmic re-
flexion studies. Manuscript, Bátatom Ltd., Budapest, BA-03—06.
Rouhiainen P. & Pöllanen J. 1999: Difference flow measurements in
boreholes KA2865A01, KA3065A02 and KTTX5 at the Aspo
HRL. International Technical Document, ITD-99-25, Swedish
Nuclear Fuel and Waste Management Company, Stockholm,
Sweden.
Sawada A., Uchida M., Shimo M., Yamamoto H., Takahara H.T. &
Doe T.W. 2000: Non-sorbing tracer migration experiments in
fractured rock at the Kamaishi Mine, Northeast Japan. Engi-
neering Geol. 56, 75—96.
Snow D.T. 1965: A parallel plate model of fractured permeable me-
dia. Ph.D. Thesis, Univ. California, Berkeley, 1—331.
Szongoth G. & Galsa A. 2003: Complex interpretation of flow and
temperature measurements. Manuscript, Bátatom Ltd., Buda-
pest, BA-03—82.
Terzaghi R. 1965: Sources of error in joint surveys. Geotechnique
15, 287—304.
Tóth Gy., Mező Gy., Benedek K. & Takács T. 2003: Hydrogeologi-
cal characterisation of blocks based on numerical modelling.
Manuscript, Bátatom Ltd., Budapest, BA-03—25.
Uchida M., Sawada A., Senba T., Miyoshi T., Shimo M., Yamamoto
H., Takahara H., Doe T.W. & Cladouhos T. 1999: Geological and
hydrological investigation and mass transport study in a fractured
system at the Kamaishi Mine. Proceedings of an International
Workshop of the Kamaishi in situ Experiments, Japan Nuclear
Cycel Development Institute, JNC TN7400 99—007, 77—86.
Van Golf-Racht T.D. 1982: Fundamentals of fractured reservoir en-
gineering. Dev. Petrol. Sci., 12, Elsevier, Amsterdam, 1—710.
Zilahi-Sebess L., Mészáros F. & Szongoth G. 2003: Characterisa-
tion of fracture zones in granite, based on well-logging data at
the Üveghuta site. Ann. Report Geol. Inst. Hung., 253—266.