Can effective porosity be used to estimate near-well protection zones in fractured chalk?

Protection of fractured carbonate aquifers is often based on a single-porosity description of a dual-porosity system. However, it is difficult to assess a trustworthy value of the effective porosity based on scientific principles; thus, a range of estimates is often suggested. The complexity of the problem is compounded by the fact that the effective porosity may be scale-dependent. This paper investigates whether it is possible to describe solute transport in fractured carbonate rocks with an equivalent porous medium model using a constant value of effective porosity. It is assumed that the dual-porosity model provides an acceptable description of transport mechanisms in fractured porous rock and that it is possible to estimate the parameters needed in the single-porosity models from results generated by the dual-porosity model. The effective porosity is estimated from the dual-porosity results that are used as targets. For Danish chalk, an effective porosity of 13% (11–17%) is estimated. However, it is demonstrated that the estimated effective porosity is only valid at the specific transport time (1 year) from which simulation results of the dual-porosity model were extracted. The effective porosity is shown to increase with travel time until equilibrium conditions are realised between the fractures and matrix, following which, the effective porosity equals the matrix porosity and will maintain this value at larger transport times. Assuming that the dual-porosity model provides a trustworthy description of solute transport in fractured chalk and limestone, a method to estimate the effective porosity of an equivalent porous medium model is presented.


Introduction
Fractured carbonate rocks form important groundwater reservoirs in many parts of the world (Medici et al. 2022), including England in the UK (Maurice et al. 2021) and Denmark (Vangkilde-Pedersen et al. 2011).The north-west European limestones form a basin that covers parts of England, France, Belgium, the Netherlands, Germany, Sweden, and Denmark.The onshore margin of this basin is today used for water supply in most of these countries (Downing et al. 1993).Chalk and limestone aquifers supply about one-third of Denmark's potable water (Vangkilde-Pedersen et al. 2011); however, fractured aquifers are confronted with a special challenge when dealing with groundwater protection because of the high flow velocities in the fractures.Two main concepts have been brought forward in groundwater protection planning and management of source water protection zones.The first one is the "capture zone" (CZ), which is defined as the entire area of an aquifer that contributes to a pumping well (Keely and Tsang 1983).The purpose of the capture zone is to protect the aquifer from contaminants that enter the groundwater system as a result of land use activity.The second one is "source protection zones" (SPZs, which are defined around large groundwater abstraction sites.The purpose of SPZs is to provide additional protection to safeguard drinking water quality by constraining the proximity of an activity that may impact a drinking-water abstraction well.Within a SPZ, there is often a distinction between zone 1 (an inner protection zone, defined by travel times of 50 days or less based principally on biological decay criteria), and zone 2 (an outer protection zone, defined by a 1-year or 400-day travel time from a point below the water table to the abstraction well).The travel time is derived from consideration of the minimum time required to provide delay, dilution and attenuation of slowly degrading pollutants (EA 2019).The focus in the following is on the SPZ outer protection zones defined by a 1-year travel time in fractured carbonate rocks, referred to as the BNBO zone (Danish abbreviation for "near-well protection area").
The area covering the 1-year travel time cannot be delineated by measurements directly (Worthington 2023).Flow and transport models are therefore required as a tool to quantify the distance corresponding to a travel time of 1 year.This task is relatively straightforward for an unconsolidated porous medium.However, for a fractured medium where fast transport is expected in the fractures and almost stagnant water is expected in the matrix, this poses a challenge.Relatively few modeling studies of solute transport in the fractured carbonate rocks of north-western Europe have been published (Bibby 1981;Brettmann et al. 1993;Little et al. 1996;Mosthaf et al. 2018).It has often been demonstrated that it is difficult to describe transport behavior in fractured porous formations using an equivalent porous medium (EPM) approach (Worthington 2015(Worthington , 2022)).To describe a fractured porous medium using an EPM model, up-scaled parameters, also referred to as "effective" parameters, should be defined.It is difficult to estimate the effective parameters that appear in the upscaled descriptions of flow and transport including effective hydraulic conductivity for flow, effective porosity and effective dispersivity for transport.Both effective dispersivity (Gelhar et al. 1992;Zheng and Bennett 2002) and effective hydraulic conductivity (Schulze-Makuch and Cherkauer 1998;Sanchez-Vila et al. 2006) are expected to be scale-dependent.The effective porosity, defined as the value of porosity of a single-porosity model that generates a correct estimate of travel time, is expected to be a timedependent parameter (Worthington et al. 2019;Worthington 2022).Additionally, to use the upscaled version of the flow and transport descriptions, it should be possible to define a REV (representative elementary volume) and the transport scale should be considerably larger than the REV (Bear 1992).It is expected that the REV of fractured chalk has to include at least several active fractures.
If the fracture-matrix system, also referred to as a dualporosity system, should be described by an EPM model that shows general validity, equilibrium between the solute concentration in fractures and matrix should be realized.Otherwise, the parameters are not constant (Worthington et al. 2019) but will be a function of travel time.This is especially the case for effective porosity, whereas the effective hydraulic conductivity is stable at smaller time scales depending on the heterogeneity structure at the site (Sanchez-Vila et al. (2006).A dual-porosity description was successful in capturing the interaction between the fracture domain controlled by advective transport and the matrix domain primarily controlled by diffusive transport (Neuman 2005).None of the transport properties, such as dispersivity or the mass transfer coefficient, were based on a hydrogeological characterization of the actual formation, including fracture distribution, spacing, connectivity, or apertures, and matrix diffusion coefficients.This is partly caused by the lack of data about these parameters, and partly because of a lack of knowledge about the relation between the fracture/matrix characteristics and the macroscopic transport parameters.Since this relationship is not established, it is difficult to assess the reliability of the calibrated continuum models; however, the authorities may require that the EPM approach is used, even for fractured aquifers.This is explained by the wish to deal with only one method for quantifying the risk of contamination.It might be more convenient to work with the same type of parameters and the same kind of models, such that the same model setups can be used for all areas and mutually compared, irrespective of hydrogeological conditions at the site.
The objective of the study is to investigate if near-well protection zones can be simulated with equivalent porous media (EPM) models where the temporal scale is a 1-year solute transport time and distance.Furthermore, the study aims to map, if data availability allows, typical values of effective (bulk) hydraulic conductivity, effective porosity, fracture and matrix porosity, and the mass transfer coefficient used to model near-well protection zones of fractured carbonate aquifers.

Theory
Modeling of flow and transport in fractured porous media may be carried out using three different approaches: The discrete fracture method (DF), the dual porosity method (DP) or the equivalent porous medium method (EPM).The three methods are described by, e.g., Bear (1992) and illustrated in Fig. 1.
In a discrete fracture model, the individual fractures and their connection are described explicitly.This approach is unattractive for two reasons-first, it is computationally demanding if the scale of interest is large compared to the fracture density, and this is typically the case in BNBO investigations; second, it is problematic to estimate the location of the fractures and the hydraulic conductivity of the fractures, as this requires detailed knowledge about the fracture location and aperture distribution.Therefore, the discrete fracture approach will not be considered in the present study.
Instead, it is assumed that flow and transport in a fractured medium can be described using a dual-porosity approach (DP, Fig. 1) as this does not rely on detailed information on the fracture configuration.Additionally, the hydraulic conductivity (permeability) of the matrix of carbonate rocks is normally low enough to ensure that the migration of solutes in the matrix is controlled by diffusion.In the fractures, advective transport will usually dominate; therefore, it is valid to describe the fractured aquifer as a dual-porosity medium, where advective transport is restricted to the fractures, and diffusion is the only transport mechanism accounted for in the matrix.
Conservative solute transport in a dual-porosity medium is described by a transport equation for each domain.The equation describing transport in the fracture zone in 2D is given by wherec f and c m are the concentrations in the fracture and the immobile matrix zone [kg/m 3 ], n f is the porosity of the mobile fracture zone [1], D ij is the hydrodynamic dispersion coefficient [m 2 /s] assuming D ij = v i α Lj + D diff where i and j designate coordinate directions.D diff is the molecular diffusion coefficient.α is the dispersivity.In the fractures, the flow velocity is normally high; therefore, advective transport is relatively more important than dispersive transport.For solute spreading in an open fracture, mechanical dispersion, given by velocity × dispersivity, will dominate over molecular diffusion, again caused by the high velocity in the fracture.v i is the flow velocity in the fractures, given by: where K f and n f are the fracture hydraulic conductivity and the fracture porosity, respectively, and J i is the hydraulic gradient.The mass transfer coefficient, β [s -1 ], controls the exchange of solute between the fracture and the matrix domains, together with the concentration gradient between the two zones.Equation (3) describes the exchange between the matrix and fracture domains which is governed by the (1) concentration difference between the domains and a (constant) exchange coefficient (Sudicky 1990): where n m is the matrix porosity [1].From Eq. ( 3) it is seen that an increase in concentration in the matrix zone is proportional to the difference in concentration between the two zones.
In the last approach, the equivalent porous medium (EPM) approach (Fig. 1), solute transport in the fractured rock is described as one united porous medium using one set of parameters.
The EPM approach requires that equilibrium conditions exist between the fractures and the matrix concerning solute concentration.Solute transport in fractured carbonates is in many cases controlled by diffusion between the fractures and the matrix.Beyond a certain transport time (at a given distance), transport into the matrix takes place at the same rate as transport from the matrix to the fractures.This is here referred to as the equilibrium transport time (t eq ).
The application of the EPM model also requires that effective (upscaled) parameters can be defined, e.g., effective hydraulic conductivity or effective porosity.Effective parameters are upscaled quantities that can capture the impact of small-scale heterogeneity on large-scale flow and transport.Hence, effective parameters describe the combined effect of heterogeneities, in the present case fractures and matrix, on the flow and transport.K eff is expected to increase with the scale of the analyzed system; however, K eff in densely fractured aquifers will reach a plateau at a relatively small scale which indicates that the hydraulic conductivity is constant at larger scales (Schulze-Makuch and Cherkauer 1998;Schulze-Makuch et al. 1999). (3) Fig. 1 Three conceptual approaches to model solute transport in fractured porous media The effective porosity (n eff ) is expected to be a function of the time scale and the transport properties of the problem at hand (Worthington et al. 2019), see Fig. 2. If the time scale is small and the exchange of solute between fracture and matrix is relatively slow, the transport through the fracture network will dominate and the fractured medium can be described as an equivalent porous medium, where the effective porosity equals the fracture porosity with an expected value less than 1% (Worthington et al. 2019).At longer time scales, molecular diffusion between fracture and matrix will be more significant and will affect the solute concentration in the fractures as well as in the matrix.After a sufficiently long transport time, it may be assumed that the diffusion from the fractures to the matrix is fully developed and the effective porosity will therefore be at its maximum, corresponding to the matrix porosity.This equals the equilibrium condition described in the preceding (t eq , Fig. 2).It should be noticed that the curves in Fig. 2 only serve as an illustrative example.Up to a given distance from a source, the EPM model with effective porosity equal to the matrix porosity can be applied at times t ≥ t eq .At time less than t eq , the effective porosity exhibits a transient behavior (Worthington et al. 2019), and it is necessary to use a dynamic description of effective porosity with values in between the fracture and the matrix porosity.This implies that the effective porosity will change when the travel distance or time, over which the transport is taking place, changes.As a result, it is problematic to use the EPM method in intermediate situations since the value of the effective porosity will change if the transport time changes or the conditions at the site change and it is not known how much it will change.
To apply the upscaled descriptions, DP or EPM, it should be possible to define a REV (representative elementary volume).At scales above the REV the effective parameters (K eff , n eff and α eff ) are constant with averaging volume while at smaller scales the parameters depend on the scale at which the parameters are defined.The model domain should be significantly larger than the REV; for a fractured medium, a REV should include several active fractures.If the spacing between fractures is, e.g., 1 m, the scale of the REV may be in the order of 10 m, and the model domain should have a size of maybe 100 m.

Methods
The study includes a review of existing information on parameters used to model solute transport in fractured chalk aquifers within the Danish chalk province that subsequently serves as input for conceptual modeling of near-well protection zones.The parameters are found using different methods: measured directly (in situ and laboratory), hydraulic and tracer tests in boreholes, and parameter values estimated from solute transport modeling.The parameters are the EPM model parameters of effective (or bulk) porosity (n eff ) and hydraulic conductivity (K eff ), whereas the DP parameters are fracture porosity (n f ), matrix porosity (n m ), mass transfer coefficient (β), and longitudinal dispersivity (α L ).Transverse dispersivity is defined as α T = α L /10.Estimated parameters and intervals are compared to values available for comparable chalk deposits (primarily from England).The outcome of the review of existing parameters (data sources and parameter values) of the relevant chalk units is included in section 'Data analysis of chalk and limestone parameters'.
Conceptual solute transport modeling is performed to analyze the impact of different flow and transport parameters on solute transport in fractured porous media within the spatial and temporal frame of a BNBO, i.e., 1 year of transport time.The conceptual solute transport modeling is carried out in three steps-(1) analyzing the sensitivity of relevant parameters for solute transport, (2) estimating transport distances for several relevant Danish chalk and limestone

Model setup
The conceptual modeling of the BNBO is carried out with MODFLOW to simulate laminar groundwater flow and MT3D to simulate nonreactive (conservative) solute transport in dual-porosity and equivalent single-porosity model conceptualizations.MODFLOW is the United States Geological Survey's modular hydrologic model.MT3D is a finite-difference groundwater-mass-transport-modeling software, often used with MODFLOW. Figure 3 illustrates the different elements of the conceptual models.The length and widths of the models vary according to the parametrization, e.g., simulations with high hydraulic conductivities require long travel distances.Different model sizes were used because larger model domains are more computationally demanding than small models with the same cell size (vertical and horizontal discretization).All models are resolved using a 10 m × 10 m grid in the horizontal plane and a constant 2-m vertical discretization.The models simulate groundwater flow in a confined aquifer (no recharge) with an 0.001 (one per mille) gradient determined by specified head boundary conditions to the left (upgradient) and right (downgradient).In the downgradient part of the model domain, an abstraction well is defined with the MODFLOW Well package with a total abstraction between 10,000 and 1,000,000 m 3 /year in the sensitivity analysis and 500,000 m 3 /year in the other simulations.The spill (simulated pollution) is always located in the upper numerical grid layer in the same row (the centre row) as the abstraction well.A homogeneous hydraulic conductivity field is assumed, e.g., the parameter values are the same for the entire model domain.
The spill is simulated as a continuous source of pollution where a solute concentration is kept constant (C 0 = 1,000 μg/L) in the spill cell throughout the entire simulation period.This was done because initial scenario runs showed that with some parameter sets a noncontinuous source will never get above the defined breakthrough concentration at the abstraction well.Also, the setup was to simulate a "worst-case spill" directly on the top of the chalk aquifer, for instance, a spill from chemical transport that accidentally discharges high concentrations of solutes in vulnerable areas of the aquifer (where there is direct contact from the surface to the aquifer).The spill is simulated as a conservative tracer, i.e., no degradation or sorption is included.Dilution only occurs as a result of dispersion and diffusion and as a result of the exchange between the fractures and matrix.Breakthrough of pollutant, the first arrival, is defined as the time when the concentration at the abstraction well reaches a value of 0.1 μg/L.Thus, the first arrival concentration (C) is defined as 1:10,000 of the spill concentration (C/C 0 = 10 -4 ).The conceptualization of this situation is affected by the local country legislation saying that near-well protection zones are defined conceptually and in practice as 1-year travel times "in the aquifer" that is being abstracted from.The concept of 1-year travel times is therefore not directly comparable with the concept of a capture zone for a well.Also, the 1-year travel distance defined in this study is only the travel time in the chalk aquifer that is being pumped.The equivalent travel distance that corresponds to 1-year travel time should be considered as a maximum travel distance because it is located directly upgradient of the well and the solute transport is conservative, so concentrations are not retarded or reduced by, e.g., sorption or degradation.The EPM and DP models have the same setup except for the matrix domain where the DP models has both fracture and matrix porosities, whereas the EPM has one bulk effective porosity for the entire domain.The DP models also included the mass transfer coefficient to describe the solute exchange between matrix and fractures.

Dual-porosity model parameters
Several processes and associated parameters affect the model-simulated migration of a solute in a fractured (dualporosity) medium.For advective transport, hydraulic conductivity controls how fast the water flows from the higher to lower hydraulic head.The hydraulic conductivity (K) in aquifers is usually obtained from pump tests.In many circumstances, the vertical component is smaller than the horizontal, for example, as a result of horizontal sediment structures, referred to as horizontal to vertical anisotropy.The horizontal hydraulic conductivity is assumed to be 3× higher than the vertical.Horizontal isotropy implies that the same hydraulic conductivity is assumed in both the longitudinal and transverse horizontal directions (x and y).Concentration at a given point is also affected by dispersion of the solute by sediment heterogeneities on different scales.The dispersion components are defined by longitudinal dispersivity (α L ) and a fixed ratio between the longitudinal, vertical and horizontal transverse dispersivities.When simulating solute transport with a dual-domain model, here described by a dual-porosity (DP) model, both a fracture and a matrix porosity must be defined.To describe the rate of this exchange in a DP model, a coefficient of exchange of the solute between the fracture and matrix domains is needed.In the MT3D model applied here, the mass transfer coefficient (β) is used for that purpose.

Assessment of effective porosity from modeling (n eff )
The first aim is to model and estimate 1-year travel distances with typical parameter values in the conceptual DP models.The second aim is to test if the estimated 1-year travel distances can be replicated using an EPM model, where the porosity is described as an effective porosity, in contrast to the DP models where separate fracture and matrix porosities are applied.
The results in terms of the 1-year travel distance and the breakthrough curve from the DP models are considered the "truth" and used to calibrate the EPM models where the effective porosity is estimated to obtain the best possible match between the two models.
According to the theoretical considerations described in the section 'Theory', the effective porosity changes continuously after the tracer has been released.At time zero, the effective porosity equals the fracture porosity, which grows continuously until a certain time (equilibrium time, t eq ) and distance after the spill, where the effective porosity equals the matrix porosity.When the focus is on near-well protection zones, defined as the 1-year travel distance, it is important to know if the estimated effective porosity reaches the matrix porosity (t eq ) within the simulated time and distance of the EPM models, or if the effective porosity is located in-between fracture and matrix porosity (t eq is not reached).If t eq is not reached, the estimated effective porosity only applies under the exact conditions, e.g., transport time and parameters, at which it was estimated.
The modeling steps to estimate the effective porosity are: 1. Set up the DP model with parameters found, e.g., from the literature review.2. In the DP model, iteratively change the location of the spill site (in the topmost grid layer) until the timing of the breakthrough (concentration < 1:10,000 of spill concentration) matches 365 days at the abstraction well in the model (Fig. 3).This is done by analysing the breakthrough curve at the well, for every simulation with different spill sites.3. The 1-year travel distance is then estimated to be the horizontal distance between the numerical grid cell of the spill site and the well, respectively.4. The spill site found by the DP model (to match a breakthrough at 365 days) is specified in the EPM model, and the effective porosity is then adjusted to match a breakthrough at the well after 365 days.
The precision of this exercise is controlled by the models' cell size of 10 m.Hence, an exact breakthrough at 365 days was not always possible to obtain in the DP model, so the spill location (numerical grid cell) that had the closest match to a responding breakthrough at the well was chosen.The effective porosity was estimated in the EPM model with two decimals of precision (corresponding to defining porosity as an integer of percentage).

Sensitivity analysis
Before the simulation of solute transport with the recommended parameter values from the literature review, a sensitivity analysis of relevant parameters for solute transport in a DP model was performed against the 1-year travel distance.The horizontal hydraulic conductivity, matrix porosity, fracture porosity, longitudinal dispersivity and the mass transfer coefficient were tested with parameter ranges found initially during the literature review (Table 1).Parameter values shown in italics were used when one parameter range was tested in the one-way sensitivity analysis.The impact of abstraction rates was also tested because advective transport increases as a result of abstraction, especially close to the abstraction well.
The sensitivity analysis was performed in the same iterative way as the estimation of the 1-year travel distance in the DP models (point 2 in modeling steps).For a given set of parameters (Table 1), a 1-year travel distance was found by simulating the spills at different locations upgradient of the well.For some of the parameter values, e.g., the high hydraulic conductivities, the models were extended in the horizontal directions (y and x, Fig. 3) to capture the entire transport path.
To compensate for the spatial discretization error caused by the 10-m horizontal grid discretization that results in a breakthrough not exactly at 365 days, the travelled distance is normalized to 365.For instance, in a situation where the horizontal distance between the spill site and the well is 210 m, and the first arrival is observed at the well on day 340, the estimated 1-year transport distance is (210 m/340 days) × 365 days/year = 224 m/year.This process is repeated for each set of parameters.

Data analysis of chalk and limestone parameters
A short introduction to the characteristics of the chalk and limestone formations is given in the following.

Geological setting/geological units within the study area
The limestone aquifers were deposited during the Upper Cretaceous (99.6-65.5 million years [Ma]) and Paleocene (Danian and Selandian, 65.5-58.7 Ma).The Upper Cretaceous marine chalk deposition in Denmark (in Danish: Skrivekridt) occurred in a part of the large North European carbonate sea and consists of at least 450 m of chalk deposits (Surlyk et al. 2013).The chalk is a fine-grained deposit where more than 80% of the chalk mass has a grain size of less than 5 μm.The chalk is a white and yellow white micrite chalk that is generally soft except for thin yellow hard grounds, occasionally black chert layers, and thin marl layers.The main component is coccoliths with small clay content (Håkanson et al. 1974).The middle Danian limestone consists of several limestone types as Bryozoan reef limestone (in Danish: Bryozokalk; Surlyk et al. 2006) and muddy or micritic limestone (in Danish: Slamkalk) here referred to as Danian chalk.A late Danian limestone (in Danish: Kalksandskalk) is deposited and classified as a calcarenite.Here, this unit is referred to as København limestone, which consists of a yellow-white limestone with some thick chert layers and is often a hard massive limestone (Thomsen 1995;Jakobsen and Klitten 1999).Overlying the Danian, early Selandian (lower Paleocene) sediments were deposited after a period of erosion and often beginning with glauconitic conglomerate and greensand of the Lellinge Greensand Formation (in Danish Grønsandskalk) followed by a massive light grey silty marl with strong bioturbation and without layering named the Kerteminde Marl (in Danish: Kerteminde mergel) (Heilmann-Clausen 1995).Figure 4 shows where the six carbonate rock types are observed in boreholes used for groundwater abstraction in Denmark.

Parameters for chalk aquifers
The existing literature on nonreactive solute transport in fractured carbonate rocks of Denmark was reviewed to investigate differences in flow parameters with rock types.The locations of available data for the different carbonate rock types in Denmark are shown in Fig. 4. The literature review focused on collecting the following parameter values from field and lab studies: effective porosity (n eff ), matrix porosity (n m ), fracture porosity (n f ), hydraulic conductivity (K), and the mass transfer coefficient (β) between the matrix domain and fractures for nonreactive solutes.The review of internationally published literature, as well as master's theses from Danish universities, showed that 34 experimental field and laboratory studies have been carried out at 19 field sites.Porosity data were reported from 20 of these studies (six studies with fracture porosity data), and mass transfer coefficients from six studies, Table S1 in the electronic supplementary material (ESM).The review shows that chalk and limestone sites from eastern Denmark dominate the collected data, whereas a very sparse data set was available from studies of porosity values at study sites in the western part of Denmark.

Hydraulic conductivity
A data set of transmissivity from ~9,500 boreholes, all obtained from pumping tests in fractured chalk and limestone wells, was collected from the national well database (Jupiter) for groundwater and water quality (Kidmose et al. 2022;Fig. 4).Transmissivity values are estimated from pumping tests carried out in boreholes that were drilled for groundwater abstraction.In general, the transmissivity data are of varying quality and were calculated from pumping tests of different pumping periods using water level measurements from either the abstraction well only or from both the abstraction well and observation well.The geographical coverage of boreholes with pumping tests is most dense in the Metropolitan area of Copenhagen and in the center of the chalk and limestone aquifers in Jutland.Much lower density occurs in the peripheral areas of the chalk and limestone groundwater aquifers in most parts of the country.Hydraulic The upper 2.5% of the estimated K values represent fracture/fissure-dominated rapid groundwater flow (~200 boreholes).Approximately 2-3% of all boreholes (~350 boreholes) have an estimated K value lower than 1 × 10 -6 m/s, which most likely represents matrix-dominated flow.The normal range for fractured limestones is defined as 10 -6 m/s < K < 10 -2 m/s with a median K value of 1.34 × 10 -4 m/s.Based on 2,100 pump tests in the Chalk aquifer in England, MacDonald and Allen (2001) found a median transmissivity of 540 m 2 /day, corresponding to a hydraulic conductivity of 1.25 × 10 -4 m/s if an aquifer thickness of 50 m is assumed.This is very close to the Danish estimate.
An analysis of the pump test data was made to determine if there is a relation between the hydraulic conductivity and aquifer test volume that the individual pump tests represent (Fig. S2 in the ESM).Data on the duration and capacity of the pump tests was available in the Jupiter database, thus the aquifer volume could be calculated for each of these pump tests.Figure S2 in the ESM shows the K value and estimated aquifer test volume of each pump test divided into matrix-dominated chalk (60 wells), 'karst' dominated chalk (44 wells), with the rest of the pump tests falling in the normal domain of fractured chalk (see Figs. S1 and S2 in the ESM).Hydraulic conductivities show a variation of eight to nine orders of magnitude.Aquifer volume varies six orders of magnitude with matrix-dominated chalk of 10 -1 -10 3 m 3 , normal fractured chalk of 10 -1 -10 4 m 3 , and the 'karst' dominated chalk of 10 1 -10 3 m 3 .The analysis indicates that there are one to two orders of magnitude difference in aquifer volume between pump tests carried out in matrix-dominated chalk and in 'karst' dominated chalk.There is no correlation between the K value and aquifer volume in the three populations, so no scale dependency can be detected in these data.In addition, the pump test data from Jupiter show much higher variability in both K and aquifer volume than reported by Schulze-Makuch and Cherkauer (1998).

Fracture and matrix porosity
The literature review primarily provided information on the matrix porosity in the relevant carbonate rock types.Only very limited information is available on fracture porosity in Danish fractured chalk and limestone.
Table S1 in the ESM shows data on porosity from København limestone, Bryozoan limestone and Upper Cretaceous chalk.The matrix porosity (n m ) varies somewhat between the three aquifer types-10-40% in København limestone; 5-45% in Bryozoan limestone; 20-50% in Upper Cretaceous chalk.Fracture porosity (n f ) is estimated to range between 0.02-1.6% in Bryozoan limestone and a single study shows a 4% fracture porosity in the Upper Cretaceous chalk.No fracture porosity data are available from Danian chalk in the reviewed literature.All fracture porosities are estimated based on the modeling of tracer test results from the experimental field sites in eastern Denmark.For comparison, according to international textbooks, the matrix porosity of carbonate aquifers in Canada, the USA, Mexico and England varies between 2.4 and 30%, while fracture porosities of 0.01-0.1% are presented (Ford and Williams 2007).Freeze and Cherry (1979) indicate that the porosity (matrix) falls in the range of 0-20% in limestone (dolomite) and 5-50% in karst limestone.The matrix porosity values of the Danish chalk and limestone aquifers shown in Table S1 in the ESM are relatively high compared to carbonate aquifer values elsewhere (Nilsson et al. 2022).

Mass transfer coefficient
The mass transfer coefficient (β) has only been estimated indirectly from modeling tracer tests.The literature review clarified that there are only limited data from Denmark on the mass transfer coefficient.Field data are available from Bryozoan limestone locations in Hellested at Stevns and in Karlstrup (both eastern Denmark) that provide β values between 2 × 10 -7 to 8 × 10 -3 days -1 and an average β value of 3.9 × 10 -3 days -1 .In the Upper Cretaceous chalk, data from three field locations are found-Sigerslev at Stevns, Marielyst at Falster, and in the Mjels chalk quarry south of Limfjorden-with information on the β values with a range of 1.7 × 10 -5 to 1 × 10 -3 days -1 and an average β value of 3.1 × 10 -3 days -1 .Because of the limited data on β with values for only the Upper Cretaceous chalk and the Bryozoan limestone, an average of 3.5 × 10 -3 days -1 is used in the conceptual modeling.The statistics of the model-derived β values based on field tracer tests in Bryozoan limestone and Upper Cretaceous chalk are given in Table S1 in the ESM.

Sensitivity of parameters
The results of the one-way sensitivity analysis clearly illustrate that the hydraulic conductivity of the aquifer has the highest impact on the 1-year travel distance.Figure 5 highlights this result with an illustration of the simulated distances (x-axis) for each tested parameter value (y-axis).The plot for hydraulic conductivity illustrates that the highest values result in travel distances above 500 m.The highest K value of 590 m/day (6.83 × 10 -3 m/s) results in a travel distance of 14,457 m.The second-highest value, K = 86 m/ day (1 × 10 -3 m/s), has a travel distance of 447 m which indicates a rapidly increasing travel distance with K values above 90 m/day.
The impact of the mass transfer coefficient on distance is complex.As described in the section 'Theory', this exchange coefficient determines the rate of exchange of solutes between fracture and matrix.The sensitivity of the parameters is high because, with high values, the solute is quickly removed from the fractures to the matrix and thereby the contaminant plume in the fractures becomes diluted, the first arrival becomes slower and the travel distance shorter.On the other hand, if the exchange between matrix and fracture is very low, the fracture concentration is not diluted and therefore a high concentration can be maintained in the fracture which leads to a longer travel distance.
The effect of matrix porosity on travel distance can be understood when considering the matrix as a reservoir for diluting the spill.If the reservoir is large, e.g., porosity of 0.4 (40%) or 0.5 (50%), it has a high potential to dilute the solute concentration in the fracture.In contrast, with a small matrix porosity, the low volume of water in the matrix has a lower potential to dilute the solute concentration of the fracture.For these processes to take place, a sufficiently large mass transfer coefficient must be specified.This is because the spill is defined as a fixed concentration in spill numerical cell.If the porosity increases, the mass of the chemical also increases, and therefore the matrix porosity affects the transport of solute concentrations.
The fracture porosity affects the travel distance, however, in two opposite directions.One process is a result of the advective transport processes where lower fracture porosity results in higher flow velocities and therefore greater travel distances.The other process reduces the travel distance when porosity decreases because the concentration of a solute in a small fracture is more easily lost to the matrix with a small fracture volume.Again, the mass transfer coefficient will affect the exchange between fracture and matrix.The combination of the two processes, working in opposite directions, results in the relatively and somewhat unexpected small effect of fracture porosity on travel distance.

Estimation of effective porosity
Based on the sensitivity analysis, it is important to quantify the impact of hydraulic conductivity on travel distances and to investigate how results from the DP model can be replicated by the EPM models with an effective porosity.The observed variation of hydraulic conductivity of the Danish carbonate rocks is described by a median of 1.34 × 10 -4 m/s (n = 8,300, see Table S2 in the ESM), a lower value of 3.60 × 10 -5 m/s (the 25th percentile) and an upper value of 5.11 × 10 -4 m/s (the 75th percentile).These values represent the different Danish carbonate rock types.Figure 6 illustrates simulated breakthroughs at the pumping well with three DP models (each having one of the K values mentioned previously).Breakthroughs simulated with the three DP models are shown as solid lines and the matching EPM models are shown as dashed lines.To match the DP model results, effective porosities of 11, 13, and 17% are needed in the EPM models with hydraulic conductivities of 3.60 × 10 -5 m/s, 1.34 × 10 -4 m/s and 5.11 × 10 -4 m/s, respectively.Effective porosities were found iteratively by trial and error as described in section "Assessment of effective porosity from modeling (n eff )".The 1-year transport distances for the 25th percentile, the median, and the 75th percentile models were 147, 184, and 321 m, respectively.Even though the breakthrough of the EPM 75th percentile and median models seems similar (Fig. 6), the travel distances are very different.Figure 7 illustrates the flow field in terms of hydraulic head for the "DP K median" model (solid black line in Fig. 6) and the simulated solute concentration at day 366 after spill in the dual-porosity model the 1-year travel time of around 18 simulation cells of 10 m in width and length can be observed.
The sensitivity analyses (Fig. 5) showed that the mass transfer coefficient (β) is the second most sensitive parameter.The remaining parameters had less impact on the simulated 1-year travel distance.Hence, further investigation was carried out on the impact of β on travel distances of the DP models and corresponding effective porosities in the EPM models.Limited data are available on β and therefore it was decided to test values one order of magnitude higher (3.5 × 10 -2 1/day) and lower (3.5 × 10 -4 1/day) than the mean value of 3.5 × 10 -3 1/day.
The DP model with β = 3.5 × 10 -2 1/day had a 1-year travel distance of 120 m (Fig. 8).For the model with the recommended β = 3.5 × 10 -3 1/day, a travel distance of 184 m was found.Based on the sensitivity analysis this is expected as the travel distance decreases with increasing mass transfer coefficient.The concentration in the fractures decreases as a result of more interaction with the matrix (dilution) and therefore the first breakthrough is observed later (or the 1-year travel times decrease).Figure 8 shows the breakthrough curves with the early breakthrough around 1 year for the DP and EPM models with different mass transfer coefficients.
In the corresponding EPM model, an effective porosity of 26% was estimated to obtain the same timing and distance as the DP model with the high mass transfer coefficient (3.5 × 10 -2 1/day).This is close to the matrix porosity of 30%.
The simulation with the low mass transfer coefficient, β = 3.5 × 10 -4 1/day, results in a very slow exchange of solutes between matrix and fractures and enables faster transport and a higher 1-year travel distance of 480 m in the mobile dual-porosity chalk domain (DP).To obtain a similar first breakthrough in the EPM model, an effective porosity between 3 and 4% has to be used resulting in the first arrival after 324-432 days (only the 3% simulation is shown in Fig. 7).

t-equilibrium (t eq )
The conceptual modeling illustrates that the effective porosity with a 1-year travel time is somewhere between the fracture and matrix porosities with estimated effective porosities Fig. 6 Simulated breakthrough (0.1 μg/L) with DP models and estimated effective porosity of the EPM models that match the timing of early breakthrough at the same distance as in the DP models.DP K 25% perc. is a dual-porosity model with a 25th percentile value for hydraulic conductivity.All DP models had the parameters: n f = 0.8%, n m = 30%, β = 3.5 × 10 -3 1/ day, α L = 1 m.All EPM and DP models were simulated with a hydraulic gradient of 0.001 and a pumping rate of 0.5 Mm 3 / year (one million cubic meters per year) To test what is referred to in the theory as t-equilibrium (t eq ), several optimizations of effective porosity in the EPM model were carried out against results at different times and distances modeled by the DP models.This was done by changing the spill time distance beyond the 1-year perspective in the DP model and changing the effective porosity in the EPM model to match the simulated breakthrough.Figure 9 and Table 2 illustrate how the effective porosities change as time and distance extend beyond 1 year.t eq occurs around 10,000 days (27 years) with the given parameterization and definition of breakthrough.After t eq , the effective porosity has reached its maximum value, which is equal to the matrix porosity of 30%.This value will remain constant although the temporalspatial scale increases further.
At temporal scales less than 3 days, the effective porosity is estimated to be 1%, which is close to the fracture porosity of 0.8%.Figure 10 shows the complete breakthrough curves for the DP and EPM models where the effective porosities are 30%.The results show that with a substantially longer   travel time than relevant in the framework of a BNBO capture zone (with a 1-year perspective), the EPM and DP models produce identical results after t eq has been reached.The estimation of t eq is investigated using the parameters K = 1.35 × 10 -4 m/s, β = 3.5 × 10 -3 1/day, n m = 30%, n f = 0.8%, and n eff = 30%.With the chosen parameterization, a relatively slow migration of the tracer is obtained.If for instance a lower hydraulic conductivity was used in the examples in the preceding, the t eq distance becomes less than 1,450 m and t eq becomes less than 10,000 days.

Discussion
The sensitivity analysis of the dual-porosity parameters revealed that hydraulic conductivity was the most sensitive parameter and the most important parameter to quantify.Hence, transmissivities from fractured carbonate formations were extracted from the Danish well database, Jupiter.The estimated hydraulic conductivity revealed surprisingly small differences between the six rock types found in Denmark, both concerning the median value and standard deviations.Even more surprising is the similarity between the British (for England) and the Danish (for Denmark) estimates, which could indicate that it might be possible to use British estimates of other parameters, e.g., the mass transfer coefficient β for Danish conditions.The estimated value for β in the Danish tests is twice as large as the value found by Cook et al. (2012) for a field test in the chalk in England.This hypothesis requires further analysis of experimental conditions in the two countries where, in particular, the time scale of the experiments is expected to be important.The analysis would probably also require that a few additional tracer tests be carried out in Denmark.By changing the effective porosity, the EPM models were set up to match the 1-year breakthroughs, simulated by the DP models (Fig. 6).The slope of the EPM models is steeper than the slope of the DP models, which could be because the DP model dilutes the breakthrough.Initially, the solute reaches a given point (in this case the well), and hereafter the solutes are exchanged to the matrix, removing some of the solutes, which happens continuously and thereby creates a less steep breakthrough slope.The same is observed in Figs. 8 and 10.Using a larger dispersivity for the EPM models would result in a more dispersed breakthrough and ensure that more similar breakthrough curves were produced.The spreading of the simulated plume in the DP model is, to a high degree, controlled by the exchange of solute between the mobile and the immobile domains and will therefore show a larger spreading than the EPM model, when the same longitudinal dispersivity is used in the two models.However, no attempts were made to improve the match between the two models as it was not expected to affect the overall conclusions.
The transport modeling aimed to parameterize the DP models with typical parameter values for the fractured carbonate aquifers, to simulate tracer concentration in the nearwell protection zone, and to estimate 1-year travel distances based on these typical conditions; and secondly, to test if the estimated 1-year travel distances can be replicated using an EPM (single porosity) model with an effective porosity.It is assumed that the results in terms of a 1-year travel distance produced by the DP models, considered as the "truth", can be used for estimation of the effective porosity of the EPM models through calibration.
Based on the results from this study, the effective porosity of an EPM model is found to be a function of travel time, if the simulated time is smaller than a certain threshold value (referred to as the equilibrium time).Similar conclusions were derived by Worthington et al. (2019) in a study of various tests in chalk aquifers in England.Effective porosity was found to be a transient phenomenon dependent upon the temporal test scale.According to Worthington et al. (2019), the estimated effective porosity is expected to be a function of the test timescale.At time zero (t = 0), the effective porosity equals the fracture porosity and it increases up to a certain time and distance after the spill, where the effective porosity equals the matrix porosity.Concerning near-well source protection zone modeling using EPM models, it is important to know whether the effective porosity reaches the matrix porosity (at t eq ) for the simulated time and distance, or the effective porosity is somewhere in between fracture and matrix porosity (n f < n eff < n m , if t eq. is not reached).If t eq. is not reached, the estimated value of effective porosity only applies under the specified conditions, including the parameters and the time scale at which it was estimated.
These results show that near-well protection zones in fractured aquifers can be modelled with an effective porosity in an EPM framework with a value between fracture and matrix porosity.The effective porosity values of 11-17% are directly related to the hydraulic conductivities under which they are found.Based on selected studies of tracer tests in carbonate rocks, Worthington et al. (2019) found effective porosities of fractured limestones and chalk of 10 -4 -10 -3 , i.e., orders of magnitudes smaller than the matrix porosity.In the same study, effective porosities at small scale (1-20 days) and large scale (30-80 years) were presented and referred to as shortterm and long-term effective porosities of 10 -4 -10 -3 and 0.38, respectively.It is clear from these results that the estimated effective porosity varies as a function of the timescale of the test.Hence, the estimated effective porosity only applies at the same time scale on which it was estimated-1 year in the case of BNBO.It is also obvious and shown (e.g., Fig. 5) that the effective porosity for a 1-year travel time is affected by the general parameterization of the model.It is found that the hydraulic conductivity and the mass transfer coefficient are the most important.Changing the hydraulic gradient would also be a condition that would change the flow regime and simulations of 1-year travel distances.The hydraulic gradient is indirectly tested in the sensitivity analysis because changing the drawdown, because of pumping rate changes, shows less effect than hydraulic conductivity and mass transfer coefficient changes (Fig. 6).The "real" hydraulic gradient is therefore not so much a result of the actual gradient between the "left" and "right" boundary condition, 0.001, but rather a result of the responding gradient from a certain pumping rate.Also, when changing the hydraulic conductivity, the 1-year travel distance changes and, similarly, so does the hydraulic gradient between the simulated spill and groundwater abstraction (pumping well).
The way sensitivity is handled in this study is by a oneway sensitivity analysis, which means that one parameter is tested with one set of the other parameters and is a simplistic approach because, in other parts of the combined parameter space, the given parameter could show different sensitivity.For instance, the mass exchange term β could behave differently if tested with different matrix and fracture porosities than used in the one-way sensitivity analysis.The observations of matrix and fracture porosities are nevertheless empirically strong and therefore the authors believe that the sensitivity analysis on β is largely trustworthy although the one-way sensitivity approach is limiting.
In groundwater protection, an inner protection zone of 50 days travel time or outer protection zones of 5 years travel time instead of 1 year is applied and delineated in some countries depending on administrative regulations.Additionally, capture zones between 25 and 100 years of travel time may be requested.Hence, at a given abstraction well, estimation of different effective porosities (if t < t eq ) is needed which may give rise to confusion.
The results presented here have implications for solute transport within the temporal and spatial scale of a wellcapture zone or catchment.Often, parts of the well-capture zone lie beyond a temporal scale of 60 years and a spatial scale of 1.5 km.On the other hand, average groundwater ages of the abstracted water are often younger than 60 years.The results of the present study show that effective porosity for such time frames should be carefully established based on additional research, since the transfer of valid effective porosity values is fully dependent on the applied 1-year travel time relevant for BNBO.

Conclusions
The effective porosity of an EPM model is demonstrated to be scale-dependent.It increases from t = 0 where it has a value equal to the fracture porosity (less than 1%) until t = t eq., beyond which n eff becomes constant and, hence, can be used at all later times.Before t eq. the effective porosity can be estimated for a specific time (or distance) but it only applies at this specific time.This implies that if results for another time (< t eq ) are requested, a different value should be estimated and used.
A dual-porosity model was set up to describe an average fractured carbonate aquifer in Denmark.The model was parameterized with effective hydraulic conductivity, fracture and matrix porosity, and the mass transfer coefficient.Hydraulic conductivity was found not to be scale dependent and the same was assumed for fracture and matrix porosity.The estimated effective porosity was demonstrated to exhibit strong time-scale dependency.The impact of dispersivity was not investigated in the present study.
An analysis of the parameters of the dual-porosity model showed that not all parameters are equally well known.In particular, information on the mass transfer coefficient is sparse, and at the same time, is an important parameter to know if accurate results should be obtained with the model.A sensitivity analysis revealed that the mass transfer coefficient is the second most important parameter and more knowledge on the mass transfer coefficient is therefore needed.It was shown that the mean hydraulic conductivity (K) of the chalk in England is relatively close to the K values in Denmark.This gives rise to speculations on the possibility that a trustworthy relation between mass transfer coefficients in Denmark and England could be established.This is an interesting subject for future analysis of effective parameters in the fractured carbonate aquifers in Denmark as well as England.

Fig. 2
Fig. 2 Sketch of how the effective porosity varies with transport time.The real relation between time and effective porosity is unknown and this example is only for illustration, not for use.Dotted line represents uncertainty intervals

Fig. 3
Fig. 3 Numerical model setup.Grid cells are 2 m high (z), 10 m long (x, flow direction), and 10 wide (y, normal to flow direction).Flow direction is from left to right.The width and length of the models vary, but the cell sizes, height (aquifer thickness), and hydraulic gradient are kept uniform.The example illustrates a 110 × 1,000-m horizontal model domain.Spills are initiated at different locations along

Fig. 4
Fig.4The six carbonate aquifers in Denmark, illustrated by the boreholes where they were observed.Median hydraulic conductivity (K, m/s), mean LogK value and standard deviation (SD) for each of the

Fig. 5
Fig. 5 Sensitivity analysis (one-way) of dual-porosity model parameters: a hydraulic conductivity, b abstraction, c dispersivity, d fracture porosity, e mass transfer coefficient, and f matrix porosity.DP models according to Table 1

Fig. 7 aFig. 8
Fig. 7 a Simulated hydraulic head and b simulated solute concentration for the conservative tracer at day 366 after spill in the dualporosity model.Hydraulic conductivity is given by the median value for chalk observations of 1.34 × 10 -4 m/s.Otherwise parameters in

Fig. 10
Fig. 10Breakthrough curves for spills at 1,450 and 2,450 m from the pumping well resulting in a travel time beyond t eq .DP and EPM models use the same parameters and model setup (K = 1.35 × 10 -4 m/s, Q = 0.5 Mm 3 /year), the DP models use β = 3.5 × 10 -3 1/day, n m = 30%, n f = 0.8% and EPM models n eff = 30%

Table 1
Parameter values tested in the one-way sensitivity analysisItalics denote values where one parameter range was tested in the one-way sensitivity analysis

Table 2
Simulated time and distances and effective porosities