Hydrogeological Models of Water Flow and Pollutant Transport in Karstic and Fractured Reservoirs

In a literature review of the recent advancements in mathematical hydrologic models applied in fractured karstic formations, we highlight the necessary improvements in the fluid dynamic equations that are commonly applied to the flow in a discrete fracture network (DFN) via channel network models. Fluid flow and pollutant transport modeling in karst aquifers should consider the simultaneous occurrence of laminar, nonlaminar, and turbulent fluxes in the fractures rather than the laminar flow by the cubic law that has been widely applied in the scientific literature. Some simulations show overestimations up to 75% of the groundwater velocity when non‐laminar flows are neglected. Moreover, further model development is needed to address the issues of tortuosity of preferential saturated fluid flow in fractures suggesting adjustments of the size of the mean aperture in DFN models. During the past decade, DFN mathematical models have been significantly developed aimed at relating the three‐dimensional structure of interconnected fractures within rocky systems to the specific fracture properties measurable on the rock outcrops with the use of reliefs, tracer/pumping tests, and geotechnical field surveys. The capabilities and limitations of previous reported hydrological models together with specific research advancements and findings in modeling equations are described herein. New software is needed for creating three‐dimensional contour maps in fractured aquifers corresponding to the outputs of particle tracking simulations. Existing software based on the equivalent continuum or multiple‐interacting continua cannot delineate the spread of pollutant migrations affected by the tortuous preferential flow pathways that occur in DFNs.

where the fluid flow and pollutant transport in groundwater are mainly addressed in the two-dimensional (2D) or three-dimensional (3D) equivalent regular network (or lattice) of fractures, which is equivalent to the actual heterogeneous permeable rock .
Recent studies have also focused on simulating fluid flow and pollutant transport in unstructured or non-regular fracture networks where discrete two dimensional fractures are randomly (Maillot et al., 2016) or stochastically generated (Hyman et al., 2016) in a 3D rock-block formation of a fractured reservoir. In fractured and porous formations, the hybrid approach, that is, separate discrete and continuum modeling approaches, is commonly applied using a revised version of the matrix diffusion and transport (MDT) model (Therrien & Sudicky, 1996). This model has been applied to 3D saturated rock blocks of the aquifer porous matrix which also includes embedded fractures. In particular, recent discrete fracture matrix (DFM) models Ţene, Bosma, et al., 2017) use separate meshes for 3D interacting continuum domains of the fractures, and of the porous rock matrix of the aquifer, to superimposed media with their own conservation equations and constitutive laws (Berre et al., 2019). The separate domains interact through the mass-exchange source and/or sink terms. The embedded and the projection embedded discrete fracture matrix (EDFM and pEDFM) modeling approaches have been recently described by Schädle et al. (2019). The non-conforming EDFM method considers a lower-dimensional fracture representation to avoid the time-consuming construction of complex meshed domains in the DFM model (Flemisch et al., 2017). Likewise, the pEDFM approach can handle also DFM with low-permeability fractures by adding non-neighboring connections between the fracture and matrix-grid cells. However, because of the complex numerical computations and input data, the scale of 3D rock-block simulations in DFM and discrete fracture network (DFN) models is frequently constrained in a cube with border lengths less than 300-500 m (Panza et al., 2016). This cube size may not be appropriate for modeling fluid flow and pollutant transport in a productive karst aquifer, the regional scale corresponding to which exceeds 1-3 km. Several previously reported 3D simulations (Fadakar-Alghalandis, 2017;Ţene, Bosma, et al., 2017) of flow and transport in DFN were applied to rescaled cubic domains of the real fractured medium using dimensionless variables and steady-state boundary conditions with fluid pressure set as one at the inlet and zero at the opposite outlet cube faces, whereas no-flow conditions are imposed on the other four faces. The rescaled form of these 3D cubes and simplified boundary conditions do not allow easy validations of DF model outputs using field measurements. In fact, the main groundwater flow direction at the field measurement scale is not one-directional (Black et al., 2017), as is the case during a pumping well test, for instance. However, DF model calibrations can be performed via numerical experiments or considering benchmark cases (Flemisch et al., 2017;Schädle et al., 2019).
Owing to the recent theoretical advancements in modeling fluid flow and transport in fractured aquifers, appropriate answers can be provided to almost all the open questions posed by Berkowitz (2002) in his review, as for instance in terms of the required link between model simulation and field observations. In particular, the use of experimental variograms based on field estimations of the investigated aquifer can now enable reproductions of stochastic 3D fracture networks that are "equivalent" (i.e., provide the same flow) to the real geometrically complex fracture systems by incorporating spatial variability of every hydrologic parameter, that is, fracture interconnections and apertures, length and density, and dip angles and azimuthal orientations (Berkowitz & Scher, 1997;de Dreuzy et al., 2002). The equivalent fracture networks obtained at the same scale of the field measurements, without missing assumptions about fracture spacing or connections and the presence of bedding planes, can provide the required link of Berkowitz (2002; …"How can anisotropy and correlation characteristics be incorporated into percolation, connectivity, scaling and stereological analyses?") between field measurements and model equations during flow and transport 3D simulations using DFN or DFM models. However, a few of the open questions posed by Berkowitz (2002) still remain unanswered; for instance, "How can we best develop "quantitative" modeling approaches that "realistically" capture DNAPL "migration" at the field scale?" This is because the theory that provides model output in the 3D map of pollutant concentrations (e.g., quantitative results), which correspond to the actual preferential pathways of the spatial spreading using the 3D DF model, is not yet fully developed in the literature. In particular, the equivalent structured or unstructured DFN fluid flow and pollutant transport can then be simulated as in a 3D network of smooth pipes (i.e., with a constant cross-section) connecting the intersecting fractures (Cacas et al., 1990;Karra et al., 2018). The residence time of the pollutant is provided, usually, by particle tracking (PT) computation methods (Cvetkovic & Frampton, 2012;Hyman et al., 2016;Shahkarami et al., 2019) that use the estimated velocities in every channel (or pipe), which are derived from solutions of the 3D flow problem to determine pollutant 3D pathways (i.e., qualitative results) and breakthrough curves of concentrations at the outlet sections. In some DFN models (Mahmoudzadeh et al., 2013), the expected breakthrough curve of the pollutant is obtained by the analytical solutions of the one-dimensional flow and pollutant transport in a 3D orientated bundle of parallel channels, which is equivalent to the studied channel network (CN), that is, the channeling method . This solution is very simple but does not give the spatial spreading of pollutant pathways. Some researchers prefer to apply the numerical method to the meshed domains (Fourno et al., 2019) in every planar fracture by using conventional continuum models and by imposing the mass continuity along every intersection of the fracture planes. Stochastic continuum models (ESC) (Hadgu et al., 2017) may also be applied. Here, discrete fractures are replaced by permeability structures into a 3D stochastic grid, and a converged solution to the field permeabilities usually requires several realizations. Frequent DFN modeling approaches prefer to mesh (or solve) the flow and transport equation in equivalent sub-networks of fractures or pipes, namely backbones (Karra et al., 2018), where the preferential flow occurs. This is because, even at a field scale (100-500 m), the DFN discretization may require a considerable (approximately 16 million) number of nodes  and several hours of computational time.
Furthermore, many studies on fracture connectivity and fracture network percolation thresholds have been carried out on regular and irregular fracture grids. This is because the fracture connectivity into the 3D rock block affects the maximum fluid flow rate percolation throughout the rock. The percolation threshold in the 3D or 2D equivalent CNs (Black et al., 2017;Hestir & Long, 1990) should be tested prior to simulating water flows in DFNs. This is because the threshold indicates whether the relationship between the observed density of the in-rock fracture connections (i.e., number of intersections per m 2 ) and that corresponding to the "equivalent" DFN representation is appropriately defined to simulate an aggregate flow over the percolation threshold (i.e., critical value). Moreover, the threshold helps users determine whether the DFN fracture lengths and count should be increased.
However, among the large number of published studies on water flow and pollutant transport in fractured aquifers, only a few have applied the layered discrete fractures (LDF) model (Kazemi & Gilman, 1993), although hundreds of millions of people worldwide receive their drinking and domestic water supplies from layered fractured aquifers (see the www.un-igrac.org map). The ESC and equivalent continuum (EC) models are widely preferred over the less approximate LDF modeling approach owing to the simplified input required during simulations, especially at large scales (>1-3 km). The LDF (Dagan, 1990;p. 1287-1289 model is a suitable simplification of the DFN model at the field and regional scales because the fluid flow and pollutant transport are addressed in a 3D set of parallel and horizontal fractures with variable apertures, and the mean aperture in every fracture plane is the main model parameter. In fact, during the horizontal flow in bedrock karst aquifers under a low head pressure gradient (usually less than 1%) the infrequent open vertical fractures work like piezometers, that is, as discontinuous narrow ditches without significant vertical flows owing to the null (or negligible) vertical component of the hydraulic head gradient (Masciopinto et al., 2008). LDF simulations provide effective solutions showing preferential flow and pollutant pathways that cannot be provided by ESC or EC models.
In this literature review, we explain how the modeling of fluid flow and transport equations can be addressed in fractured karst aquifers, considering validation methods and quantitative uncertainty estimations of simulation results. We focus on the flow and transport mathematical model equations and the importance of the tortuosity of the preferential flow pathways on the fluid flow velocity estimation in every 2D rough fracture and each pipe of the 3D network. Thereafter, we highlight the efficacy of geophysical methods to support the DFM, DFN, or LDF model during simulations, as forecasted by Berkowitz (2002), providing data for input and model validations. Finally, some research findings suitable to upgrade fluid flow equations in current models, taking into account nonlaminar or turbulent flow, and to "realistically" capture pollutant migration in 3D contour maps are discussed.

Equivalent Fracture Pattern Generation
A variety of software for fracture pattern generation has been developed to reproduce realistic fractures of rock blocks. The available code types include: (a) random fracture generators; (b) pure stochastic fracture generators (SFGs); (c) geostatistical fracture generators (GFGs), which are not yet fully developed to handle complex inputs and time-consuming computational procedures to reproduce 3D experimental variograms and spatial fracture data interpolations; and (d) process-based generators, which account for the mechanical (i.e., rock fracturing) and hydrogeological processes (Lei et al., 2017).
However, the stochastic method often fails to capture spatial variability and connectivity of actual fracture networks due to field data limitations. Berkowitz (2002) highlighted the need to improve stochastically generated fracture patterns to incorporate specific geometric, hydraulic, and/or chemical properties of the fracture formations that can be derived from site-specific field investigations at small and large scales, as proposed in the GFG or SFG. Some researchers (e.g., Fumagalli et al., 2019) proposed conditioned mesh generation under conforming, non-conforming, and non-matching conditions at fracture intersections for preserving real fracture connections of the fracture networks where the field data are available. Although the connectivity and flow or mechanical properties of fractures are key factors in several applications, only a small fraction of fracture intersections of a whole rock formation can generally be directly observed and studied during field investigations and reliefs. Therefore, because purely deterministic fracture network generation is frequently not possible (Chilés, 2008), the stochastic and geostatistical methods are frequently applied in DF models, since GFG and SFG require minor field data input compared to deterministic or hybrid (stochastic and deterministic) methods. Moreover, using experimental variograms the CFG and SFG can indirectly incorporate the observed fractured rock features and properties, such as fracture geometry, aperture, and connections. This can be accomplished by employing the field investigation results that can be considered key elements for improving the representative (or equivalent) stochastic fracture-network generations. In particular, GFG is based on the study of the spatial covariance of the investigated fracture parameters (mainly the fracture aperture, density, length, and orientation) using experimental variograms and it is suitable for yielding representative fracture patterns that are equivalent to the studied rock formation. It uses fitted variogram models for each variable (Gelhar, 1993). SFG requires probability distributions (Darcel et al., 2009) applied to the field measurements, specifically the truncated power law distribution to determine fracture sizes, Poisson or uniform distribution for fractures spatial positions, and Fisher distribution for orientations. GFG are suitable in deterministic DF model simulations whereas during stochastic DF model simulations that are based on ensemble averages of DFN domains, multiple SFG realizations are required (Berkowitz, 2002). Thus, several successive Monte Carlo runs accounting for the conductive or non-conductive probability trends of grid node neighboring fractures (or channels) are required until the experimental values, that is, fracture length and conductance (Dessirier et al., 2018) are met. Sandve et al. (2014) proposed DFN improvements using "preconditioners", that is, the mass conservation property and real fracture geometry, for preserving the representative flow characteristics of rock formations during the simulation upscaling.
Input data for GFG and SFG could be derived from the elaboration of the fracture outcrop reliefs using FracPaQ (Healy et al., 2017) or similar codes, and from the results of field tracer/pumping tests in boreholes. Some software, for instance ADFNE 1.5 (Fadakar-Alghalandis, 2017) (Figure 1), can implement both the stochastic fracture generation and variogram methods for fracture permeability generations in a 3D DFN equivalent to a block of karst aquifer.
Hence, the DFN given by GFG or SFG can be considered an equivalent network (Huang et al., 2017) of the interconnected smoothed fractures of different apertures. However, both SFG and GFG should be applied only when a minimum of field measurements are available to obtain a realistic flow model output and reduce the uncertainty of the results to acceptable values. This is because field investigations and observations of exposed rock walls, tunnels, or ditch measurements by means of borehole tests cannot always be exhaustive of the deep groundwater flow, especially when the flow and transport properties must be scaled up to the entire regional domain of the fractured aquifer.
The fracture aperture variability in each fracture plane, or among fractures in a DFN, can be defined by a geostatistical aperture generator (GAG) based on the experimental directional variograms of the aperture covariance derived from field estimations. The mean aperture at each position or grid node of an LDF model can be computed, for instance, using data from steady pumping (or injection) borehole tests, by inverting the analytical solution of the radial flow to the well (i.e., the Thiem equation) (Masciopinto et al., 2008;Tsang, 1992) into the parallel set of LDF with a mean aperture 2b to obtain where Q w (L 3 /t) is the pumping (or injection) well flow rate, Δϕ (L) is the drawdown (or mound), r i (L) is the pumping (or injection) influence radius, r w (L) is the borehole radius, n is the effective (i.e., secondary) rock formation porosity, B is the saturated flow thickness, and μ/γ (L⋅t) is the water viscosity/specific weight ratio. A minimum number of aperture estimations (30-40) (Gelhar, 1993) should ensure a good fit of the experimental variogram covering a wide range of spatial distances at regional scales.

Guidance in the Selection of the Model
The selection of the appropriate model equations that should be applied in fractured bedrock aquifers is not easy. In this section, the modeling approaches that have been implemented in commonly applied software in the literature (Figure 2), considering their concepts, and peculiarities, are discussed comprehensively.
The geometry of embedded fractures and their sizes must be known a priori. The application of the DFN (or DFM) requires input data for fracture geometry and distributions into the rock formation, for example, dip angle and direction, azimuthal orientation, fracture size, spacing, and aperture values. This is not the case for pure stochastic models, which are based on the equivalent continuum modeling approach or ESC.

Rock Aquifers With Impermeable Rock Matrices
Fluid flow and pollutant transport throughout a saturated 3D DFN embedded in an impermeable rock matrix can be addressed using the channeling theory, which provides solutions in every single channel or stream tube of the 3D CN defined by fluid flow pathways (i.e., backbone). These can be determined in 3D DFN using the results of fracture intersections and network connectivity fields and employing the graph connectivity method, as in dfnWorks (Aldrich et al., 2017;Hyman et al., 2015;Karra et al., 2018). In the fully discrete modeling approach, after the fracture network generation, the rigorous extraction of nodes and edges of the backbone skeleton (topology) ( Figure 3) associated with the length, aperture, and any other conductivity parameters of fractures can be performed. The graph structure provides significant information, and in particular neighborhood connections, at every node. Here, the computational code for solving the 3D fluid flow problem can be applied to the pipe network, thereby ensuring mass continuity at each intersection node of the pipes (or channels) in accordance with the below equation 10.1029/2021WR029969 5 of 18 Here, the number of all inflow and outflow grid channels M at the intersection node can be determined by considering the 3D (or 2D, in a single fracture) pipe network topology. Moreover, Q ij denotes the pipe flowrate between the extremity nodes i and j, and it is usually defined for a cross-section 1 × 2b in accordance with the cubic law given by −γ (2b) 3 (12μ) −1 ∇ϕ. Equation 2 imposes the incompressible fluid-flow mass-conservation in the absence of sources and/or sink. When applied to all intersection nodes, this equation yields a system of linear equations suitable for determining the hydraulic head at each CN node. In these CN-based DFN models, the permeability inputs corresponding to the fracture apertures, fracture count, and interconnections can be generated geostatistically.
The pollutant transport in a bundle of 3D (not intersecting) channels orientated in the main flow direction was initially presented by Tsang et al. (1988Tsang et al. ( , 1996 who suggested the particle-following tracking numeral technique. This Lagrangian method provides breakthrough histograms of the collected number of particles at assigned time intervals into the observation sections (or well), as well as particle pathways of the pollutant migration into the 3D network of channels (Moreno & Neretnieks, 1993). The results provide the elapsed time required by every numerical (indivisible) pollutant particle injected in the inlet position to reach the outlet positions in a 3D CN by following any specific pathways. Specifically, particles entering a channel intersection (node) are distributed into the outlet branches following the values of the density probability function f(x), where x(t) is the vector of the spatial particle position. Hence, the particle-following tracking technique defines the pollutant pathways, that is, the spatial spreading of the pollutant owing to the advection and dispersion phenomena (Moreno et al., 1985). The pollutant concentration resulting in a target position (i.e., outlet section or a well) caused by the continuous pollutant injection can be defined as: where Δm k f (x k ) is the expected mass of the injected pollutant M 0 , which is transported by a generic particle cluster with travel time τ k provided by the PT result; nc is the number of clusters of particles having different travel times at the outlet position (i.e., histogram); and  (t −1 ) accounts MASCIOPINTO ET AL.
10.1029/2021WR029969 6 of 18  for the biodegradation (or inactivation) first order kinetic of the pollutant or pathogen removal. If a pollutant adsorbs and diffuses into the rock matrix, C 0 (ML −3 ), which is the maximum concentration that can reach the position x(t), will reduce. Specifically, using the above Lagrangian approach in the case of pathogen transport in single fractures,  can be estimated from a specific equation derived from the analytical solution of the coupled differential equations of the hydrodynamic advection-dispersion of the pollutant transport in a single fracture with variable apertures and the kinetic equation accounting for pathogen removal by physical retention or straining, inactivation, and attachment. Mathematical details are reported by Masciopinto et al. (2019;p. 569-570).
In a Lagrangian framework, the analytical PT solution of the total (passive) advective solute mass flux (or concentration), C(t, z) was obtained in the integral form Hyman et al. (2016) by where t' = (t−t 0 ), t 0 is the injection time, and H (t) is the Heaviside step function with H (t') = 1 for t' > 0, and H (t') = 0 otherwise. Equation 3a gives the cumulative distribution function (i.e., histogram) of solute passing through a control section (outlet plane) at distance z away from the inlet solute mass section at a time t. The expression 1−C (t, z), that is, the complementary cumulative distribution function, was used by the same authors to estimate the power law scaling of the tail of the solute mass breakthrough curve. Equation 3a is applied to DFNs with a large number of fractures (i.e., 200) where the pollutant dispersion owing to pathway intersections can be neglected (Cvetkovic & Frampton, 2012). To this end, Equation 3a can be applied to graph methods, thereby extracting the backbone subnetworks from DFN representations.
PT histograms of the arrival times of collected particles can also be fitted by the analytical solutions provided by the channeling theory, such as in the 3D channeling model (Moreno & Neretnieks, 1993;Tsang et al., 1988).
During particle-following tracking simulations it is assumed that there is a complete (instantaneous) and perfect mixing of the pollutant at every intersection of channels and that the pollutant-mass flow rate coming to a node is divided into the outlet branches, following probability conditioning (e.g., weighting) of the outflowing fluid flowrates (i.e., the majority of pollutant particles follow the sequence of channels with the highest velocities). Alternative pollutant mixing options according to the fluid velocity or the Peclet number at each intersection of channels were also investigated (Berkowitz et al., 1994). Kang et al. (2015) noted that a high heterogeneity of the fracture conductivity (σ 2 ln K = 5) leads to negligible effects of the streamline routing mixing option (i.e., all pollutant particles follow the adjacent pathway outflowing from the node) on the pollutant spatial spreading.
It is noteworthy that t to determine the fluid velocities along particle pathways the above-mentioned PT transport simulations must obtain preliminary solutions to the fluid-flow equations. Contrarily, the DFN meshed models (Fumagalli et al., 2019), first study the 2D fluid flow and pollutant transport for each fracture to simulate the flow and transport across all the interconnected planes constituting the 3D fracture pattern.
Thereafter, the pollutant transport in the DFN can be simulated either using PT (the Lagrangian method) or by solving the hydrodynamic advection-dispersion equation in every interconnected fracture plane of the 3D DFN (Eulerian method) using meshed domains. These model computations are simplified in a layered fractured aquifer using the LDF modeling approach. Here, the flow and pollutant transport equations can be solved in every single 2D channel network of a 3D set that is equivalent to the real layered aquifer using 2D particle-following tracking.

Fractured Aquifers With Permeable Rock Matrices
The flow modeling approach into DF embedded in a porous rock matrix was initially proposed by Sudicky and McLaurin (1992) and was defined as the MDT modeling approach. The MDT model has been applied with few but important modifications in modern modeling approaches, that is, the EDFM and the pEDFM , which account for additional non-neighboring connections between fracture and matrix grid cell. These modeling approaches consider two separate interacting media and the corresponding governing equations representing flows into the fracture cells and porous matrix domain. Here, following a hierarchical approach (Lee et al., 2001), the permeability contributions from short fractures within a grid mesh size are used to improve the average permeability of the porous matrix. The long fractures are then discretized on a separate numerical grid. Thus, porous domains of pEDFM also consider single fractures that are embedded in the porous rock matrix. The fractures overlap with the matrix cells; this ensures each porous-rock cell contains exactly one fracture cell per fracture. The fracture intersections can be used to determine the flow-transmissibility coefficients between the matrix and fracture cells as well as between interconnected fracture cells. Flemisch et al. (2017) compared the efficiencies of flow simulations performed using different discretization methods for the DFM domain. Two distinct equations are required to perform the transport simulation using DFM models. These equations represent the advective-dispersive pollutant transport in the fractured porous medium. One equation holds for the fractures while the other holds for the porous matrix. The pollutant-mass exchange between the porous and fractured domains is similar to the seepage-type pollutant-discharge exchange (i.e., source and/or sink-type) between two leaking aquifers placed adjacent to each other. Therefore, exchange in the water and pollutant mass flow rates are considered in the governing fluid flow and pollutant transport equations (Therrien & Sudicky, 1996). Additionally, the diffusive Fickian pollutant-mass exchanges between the porous matrix and the fracture fluxes are considered in the DFM models. Section 4.2 describes a detailed advective-dispersive transport equation of nonconservative pollutants into a single fracture. Furthermore, Fumagalli and Scotti (2013) have developed advective-diffusive transport equations and source functions in the EDFM coupled matrix and fractured domains, whereas Odsaeter et al. (2019) have recently presented the advective pollutant-transport equations and their solutions considering a modified EDFM model. The integration of the separated governing flow equations pertaining to the fracture domain and matrix zone must respect the flowrate (or pressure head) and pollutant concentrations imposed on the boundary, thereby ensuring mass-balance control (i.e., continuity) at each time step in the rock volume under consideration.
Because the data input requirement remains a challenging task for solving fracture-porous flow and transport, especially at the regional scale (>1 km) where there are hundreds or more fractures (e.g., about 10 9 fractures in the whole unsaturated zone at Yucca Mountain; Doughty, 1999;p. 77), in these aquifers, the equivalent stochastic generation of permeability and ESC, or multi-interacting ESC, may replace DFM modeling approaches (Reeves et al., 2008). However, in ESC model simulations, the conversion of discrete fractures or more common fracture zones to permeability structures on stochastic grids may lead to a high uncertainty of simulation results (Weijermars & Khanal, 2019) when the number of permeability measurements (i.e., input data) is low.

Hybrid Discrete → Continuous Modeling Approach
At the regional scale, many coupled DFN/EC multiscale modeling simulations in fractured rocks are commonly applied to overcome the above-mentioned scarcity of data input. Usually, to determine rock block properties at a specific site, the discrete DFN modeling approach is preferred at the locale scale (<500 m) by using hydraulic conductivities from field investigations to determine the fracture sizes. Subsequently, EC or MINC realization can be used for upscaling (Faybishenko et al., 2015) the simulation results, that is, the pollutant spatial spread into the 3D rock outcrops, to the full scale of 1-3 km to fit the observed target field concentrations.
One commonly used upscaling code (www.woodplc.com) is performed utilizing empirical formulas that provide "equivalent porosity" and "transmissivity" values (Joyce et al., 2010;R-9-20). Similarly, other software (www.golder.com/fracman) can combine detailed 3D flow solutions of DFN with the scaled-up 3D solutions given by the MINC approach.
However, it is noted that following this hybrid multiscale approach, 3D maps of the spatial pollutant concentrations in the fractured aquifers are obtained using sequential EC model simulations by stressing the mass continuity. This means that the final simulation outputs provide non-realistic concentrations (Hadgu et al., 2017) for reservoirs with few fractures and an impermeable rock matrix.

Implementation of New Flow and Transport Equations
All DF modeling approaches in the present work usually apply the cubic law to study the laminar fluid flow in fractures under steady-state conditions. Subsequent technical improvements of the transport model solutions applied to the DFN discrete models of karst aquifers have been proposed in the literature. Clemens et al. (1996), Liedl et al. (2003), and Shoemaker et al. (2008) proposed the replacement of the cubic law, which is usually applied to estimate the water flowrate using hybrid models that couple quadratic pipe flow equations, such as the Darcy-Weisbach equation (Bear, 1972;p. 126), to MODFLOW. In fact, Cvetkovic and Frampton (2012) found that the cubic law could lead to an early breakthrough curve of pollutants owing to velocity overestimations up to four times faster than those obtained using the quadratic law. The Darcy-Weisbach equation may account for laminar, non-laminar, or turbulent fluxes using the appropriate relationship between the Fanning friction factor (f) and the Reynolds number (Re); this relationship facilitates the application of the Darcy−Weisbach flow equation under any Reynolds number that may occur simultaneously in different fractures of the DFN, by where f (Re) [-] is calculated from the Reynolds number Re = 4R h U/ν, where ν (L 2 /t), U ij , and R h ≅ b denote the kinematic viscosity, is the water velocity, and hydraulic radius estimated for fluid flow in the x ij -direction into a large (regarding the aperture) rectangular cross-section of a single fracture with aperture 2b, respectively. Moreover, ϕ denotes the pressure head (L), and g (L/t 2 ) is the acceleration due to gravity. In particular, to facilitate the application of Equation 4, a specific experimental f-Re relationship is required to determine the f value in every fracture fluid flow of the DFN. For this aim, results of pumping-well tests (Masciopinto et al., 2010) at high flowrates (up to 116.5 l/s) in a fractured and karstified aquifer, and from laboratory experimental investigations carried out by other researchers (Ji et al., 2008)  For practical applications of Equation 4, hydraulic heads and water velocities given by the system of Equation 2 using the cubic law can be used in Equations 4a and 4 to determine the Reynolds number and effective velocity mathematical statement, respectively, in every fracture of the DFN. Then, solving the system of Equation 2 using the revised flowrates, the correct hydraulic head in every node, that is, accounting for losses owing to drag forces due to the non-laminar or turbulent flows in every fracture, can be obtained. Thus, correct velocities can be obtained using Equation 4. The replacement of the cubic law with Equation 4 avoids overestimations of the groundwater velocities up to 75% during the above-mentioned pumping-well tests.

Tortuosity Impact on the Flow Velocity in Karst Aquifers
MODFLOW−2005/CFP (https://pubs.usgs.gov/tm/tm6a24/) implemented tortuosity (Clemens et al., 1996) in the Darcy-Weisbach equation to study water flow in a karst conduit embedded in a fractured rock matrix of the aquifer (Medici et al., 2019).
Carman (Bear, 1972;p. 110) proposed the introduction of the tortuosity factor (<1), that is, to study fluid flow along a pathway in a pipe of non-circular shape of a porous medium. Equation 5 provides the corrective dimensionless coefficient of the flow velocity between two cross-sections with distance L, concerning the effective path length, L e . Cacas et al. (1990;p. 494) used the tortuosity factor to reduce the time residence (or to increase the velocity) of the tracer transport during the particle-following simulation in the DFN. Tsang (1984) investigated the effect of flow path tortuosity on the breakthrough curves of a tracer transported in fractures. They noted that when the contact area between the wall surfaces in the fracture was high (>35%), the mean fluid velocity and rock permeability estimated using the arithmetic mean fracture aperture 2b c in the fracture plane as derived from the cubic law can be three or more orders of magnitude less than the effective fracture velocity and permeability (Lang et al., 2014) along the actual preferential flow pathways (stream tubes). This preferential flow produces a pronounced skewness of typical non-Fickian breakthrough curves (Peclet number < 5) because of the macro-dispersion of the tracer transported in rough fractures, which is known as the "channeling effect". Tsang et al. (1988;p. 2048) defined an amplified mean "tracer" aperture 2b t to adjust the residence time of simulated breakthrough curves during simulations. Simulation results obtained in subsequent studies revealed that adjustments in the size of the mean aperture, ranging from 1.3 to 2.2 times, were required in DFN models to account for the tortuosity of the preferential tracer pathway at a low Peclet (<5) number. Minimum corrections are instead required for higher (>20) Peclet numbers.
In a later investigation (Masciopinto & Palmiotta, 2012), the tortuosity factor for water flow preferential pathways in every fracture of the LDF model was defined as: Empirical τ values can be derived from water velocity measurements given by specific tracer tests in different worldwide fractured aquifers.

Pollutant Transport Into a Single Fracture of the DF Model
More recent improvements in DF model equations (Shahkarami et al., 2019) include the non-equilibrium pollutant transport in fractures by coupling the kinetic equation of the pollutant mass adsorption-desorption and diffusion into the rock matrix or in the lateral dead-ends of fractures, to the advection-dispersion transport equation. Furthermore, for the reactive pollutant transport in fractures, a simple and efficient method was presented by Larsson et al. (2013) to account for the effects of the specific flow-wetted surface (sFWS) during pollutant transport with reactions. This sFWS is the fraction of the fracture wall contact area where pollutant adsorption occurs during the transport in fractures of DFN. In every fracture of assigned aperture 2b ij of the discrete model, the non-equilibrium transport equation of sorbing pollutants with mass-transfer during flow can be set as Here, it is assumed that the concentration of the pollutant does not modify the fluid density. The pollutant mass-transfer between fracture flow and stagnant water in the rock matrix or dead-ends of channels that are embedded in a rock matrix surrounding the flow is represented by two diffusive-type terms in the mass conservation Equation 8, that is, TD (-), the time delay coefficient, and λ d (t −1 ). The latter is the scale-invar-

Model Calibration/Validation
Many tracer experiments in boreholes that were drilled following specific dip angles into fractured rock outcrops in Sweden were described by Tsang et al. (2015Tsang et al. ( , 1996 to characterize the 3D fracture network equivalent to the investigated rock. Pumping wells and field tracer tests are usually carried out to calibrate the flow and transport of LDF models in layered karstic aquifers (Masciopinto et al., 2008). Results of tracer tests carried out in boreholes ( Figure 4) are suitable for defining the average hydraulic aperture 2b (and permeability) of single fractures in DFN at several spatial positions of the investigated fractured aquifer. This can be accomplished by inverting the flow-velocity equation, as derived by the cubic law for fractures of cross-section 2b × 1, that is, for the assigned head gradient of the horizontal groundwater flow and aquifer secondary porosity n. To this end, the specific groundwater discharge, that is, U × n, as defined by the Darcy law, can be estimated using tracer breakthrough concentrations curves (Masciopinto, 2021) into the well, where n = 1. The results of the pumping and tracer tests given by Equations 1 and 9, respectively, provide inputs for GAGs to obtain suitable fracture apertures in the DFN, DFM, and/or LDF models.
Furthermore, non-invasive geophysical investigations have been shown to be valid tools for calibrating DFN or DFM model coefficients. In the context of DFN model simulations, time-lapse geoelectric resistivity (ERT) can detect differences in rock porosity (Day-Lewis et al., 2003) and water saturation degree, which are suitable for defining the mean fracture aperture (i.e., 2b × B = n) of an aquifer thickness, by measuring anomalies in rock volumes in terms of electric resistivity changes. For instance, Figure 5 shows the agreement between the suction head contours determined by the solution of a DF model, and the rock water content obtained from ERT measurements during the unsaturated flow into a fractured rock outcrop. The wetting front shape and depth, at an assigned time, have been used to compare the model output and the apparent resistivity (i.e., water saturation) changes derived from the subsoil tomography. This comparison also suggests an effective correspondence between the modeled positions of the large fracture apertures (or cavities) and the actual positions given by ERT where the rock water content is high. Thus, the ERT MASCIOPINTO ET AL.
10.1029/2021WR029969 11 of 18 can be a suitable tool to validate DF model outputs because of the high correlation between the rock water content and electric current. In fractured rock formations, resistivity methods are effective for detecting the high variability of the fractured systems, which is caused by changes in rock electrical resistivity R ρ (Ω m) because of cavity, seawater encroachments, and clay or compact rock layer inclusions. Along the coast, ERT is widely used to visualize the sea encroachment in aquifers because R ρ decreases abruptly, reaching values less than 5 Ω m when salinized water is present in the fractures. Beaujean et al. (2014) used ERT to constrain seawater intrusion models because of its high sensitivity to salt dissolved in groundwater and its high spatial coverage. In particular, a sequential approach was used to identify hydraulic conductivity and dispersivity in density-dependent flow and transport models based on surface ERT-derived mass concentration. The azimuthal geo-electrical resistivity measurements can also be performed to detect the secondary blockrock porosity and orientations of preferential fracture flows (Lane et al., 1995;Taylor & Fleming, 1988). In this method, resistivity measurements are collected by employing a square array configuration with four electrodes placed at the corners of a square measurement area.
Moreover, the microgravity anomalies can be used to detect cavities in the karst aquifers and subsoils. The electrical and electromagnetic methods include a wide range of techniques based on the electrical and/or magnetic rock properties (Bechtel et al., 2007).
Finally, at a smaller scale, ground-penetrating radar (GPR) has also been demonstrated to be a powerful tool for obtaining data input for modeling fluid flow and pollutant transport in fractured formations. The use of GPR in detecting fractures in granitic rocks has been widely applied in the literature (Dorn et al., 2012;Tsoflias & Becker, 2008). Here, fracture apertures can be detected by resolving temporal changes of the tracer concentration throughout fractured aquifer zones.

Uncertainty in DF Model Results
The main uncertainty of model output is related to the assumptions made by investigators in the assignment of coefficients (i.e., transmissivity) of DF flow models, where fracture networks are established and their properties are obtained using GFG and GAG or using stochastic simulations. The uncertainty is caused by the absence of ergodicity of the considered stochastic variables, such as the aquifer transmissivity or MASCIOPINTO ET AL.
10.1029/2021WR029969 12 of 18 the mean fracture aperture. This is because the stochastic method requires an ergodic assumption, which ensures that the spatial values of variable covariance due to the fracture heterogeneity remain uncorrelated (i.e., stationaries) over an assigned distance (i.e., the integral scale) (Gelhar, 1993;p. 232-233). The lack of ergodicity can produce an uncertainty of predictions of DF model simulations that would be greater than the uncertainty related to the minimum data requirement (Dagan, 1990). Following pure stochastic model simulations, the model coefficients explicitly account for the natural heterogeneity of the fractured medium, using the averaged-ensemble of the DFN domain. This indicates that the averaged-ensemble model results may be derived after a large number of SFG realizations of assigned rock fracture permeabilities and lengths, for instance, 500 in the study by Dessirier et al. (2018). However, equating predicted ensemble-averages of spatial concentration moments to the experimental pollutant plume (i.e., the ergodic assumption) is possible only when the initial volume size occupied by the pollutant plume (i.e., field data) is large enough to capture all the relevant spatial correlation scales of heterogeneities that have been involved during groundwater flow and transport simulations in the studied aquifer (Sposito, 1997).
The single-realization approach, using GFG and GAG in the DF models relies on the "implicit" ergodic assumption because the resulting fracture patterns (or apertures of individual fractures) are based on the bestfit of the experimental semi-variogram derived from field measurements. In this case, every fitted semi-variogram model of the spatial covariance can be considered "equivalent" to that of the experimental variogram when the number of measurements in the latter, that is, the number of field tracer and pumping well tests, is large enough to ensure the constancy of the spatial covariance over the integral scale.
Many studies (Fiori & Bellin, 1999) have analytically estimated the uncertainty affecting the second moment (i.e., the variance) of the spatial pollutant concentrations based on the outputs of transport models. In particular, Fiori and Jankovic (2005) defined the uncertainty in the estimation of the spatial moment of pollutant concentrations when the pollutant source had a size in depth (z−direction) greater than in the x and y dimensions, which are the prevalent groundwater flow component directions, as in the case of a contaminated well. In a layered fractured aquifer, the uncertainty in the LDF model estimations is defined by the variance of the fluctuations F 11 around the mean originated by water flow velocity, and it can be set in the x−direction as: where L 1 and L 3 are the horizontal (x and y) sizes of the polluted area at the source; ℓ 1,3 is the correlation length in the x−y direction and defined by the semi-variogram of the particle trajectory x ij in the fracture plane; and X 1,1 (L 2 ) is the estimated variance of the displacement of one pollutant particle, including advection and dispersion. Thus, the non-ergodicity uncertainty can be estimated from the square root of the ratio between F 1,1 (plus the initial moment at the source) and the corresponding maximum non-ergodic variance value. If the size of the pollutant source is negligible, the maximum uncertainty is achieved using the square root of the ratio of F 1,1 (plus the initial moment value) and the corresponding maximum non-ergodic vari- (Dagan, 1990;p. 1288), where σ k 2 is the variance of the hydraulic conductivities in the case of an LDF model with a mean conductivity value K A (L/t), and U is the groundwater velocity.
Alternatively, the same uncertainty can be derived from the mean dispersion coefficient  1,1 D , which is estimated by the predicted breakthrough curve of the pollutant at the target well placed at an assigned distance L from the pollution source. This estimation can be made using PT results (Moreno et al., 1985) to be where t 0.1 and t 09 are time percentiles obtained from the pollutant breakthrough curve. Therefore, the model uncertainty of an LDF model output can be determined by the coefficient of variation of the mean dispersion coefficient (Equation 11), considering that the standard deviation of   (Dagan, 1990;p. 1285). The ergodic uncertainty was about 15-20% in largescale simulations (Masciopinto et al., 2019). Additionally, further model uncertainties owing to limited amount of available data and errors in the estimation of model coefficients (Hyman et al., 2016) must be considered to have total DF model uncertainty.

Conclusions
This review highlights improved equations with wider implications in the areas of fluid dynamic software and hydrogeological models that can be used to further study fractured karst aquifers, including new DFN, DFM, and LDF codes. These equations account for the laminar and non-laminar flow in tortuous paths in DFN as well as for the non-equilibrium mobile−immobile pollutant transport caused by mass transfer in the stagnant fluid of rock matrices or channel dead-ends.
We show that DFN and LDF models, as well as (p)EDFM models, are powerful tools to solve fracture water flow and pollutant transport in fractured and karstic rock formations. However, their full implementation at a regional scale has not yet been achieved with modern software. Currently, hybrid codes apply EC or the MINC approaches to predict maps of concentrations of the spatial spreading of pollutants released into an aquifer, although these models cannot be applied when there are few fractures. Furthermore, DFN models can provide impressive results using 3D PT simulations of the pollutant transport according to the tortuous preferential flow pathways, which occur in karst fractures and are unpredictable by EC and MINC models. However, there is not yet a specific method to convert the spatial PT outputs in 3D maps of pollutant concentrations.
Additionally, the collection of the required data input for current full-scale DFN and DFM model simulations in fractured aquifers remains very complex because of the many variables and rock properties that must be provided in the input of DF generators, and the complex computational procedures. For this, field geophysical investigations could integrate field measurements providing more input data to define the "equivalent" DF network by minimizing the uncertainties of the simulation results. The potential of geophysical field surveys has not yet been fully exploited in flow and transport simulations in fractured reservoirs. Such techniques have been recognized to be highly efficient because they provide key information, such as rock permeability, porosity, fracture orientation, and direction, as well as outputs (i.e., subsoil tomographies, profiles, or images) that are suitable for the calibration and validation of DF models. B (L) saturated flow thickness 2b, 2b t , 2b c (L) mean fracture aperture C, C T , C 0 , C s (M/L 3 ) dissolved concentration of the considered chemical D (L 2 /t) tensor of hydrodynamic dispersion coefficient of the pollutan  1,1 D (L 2 /t) mean stochastic estimation of dispersion coefficient in fractures f (-) Fanning friction factor f (x) probability function defining the particle trajectory (x) F 11 (L 4 ) variance of the fluctuations around the mean of one pollutant particle displacement originated by water flow velocity g (L 2 /t) acceleration of gravity H (-) Heaviside step function k (L 2 or darcy) single fracture permeability K (L/t) single fracture hydraulic conductivity L 1 and L 3 (L)

List of symbols
horizontal (x and y) sizes of the polluted area L (L) distance between two flow cross-sections in the fracture L e (L) effective flow path length in the fracture ℓ 1,3 (L) correlation length in the x-y direction M (-) number of converging channels in a generic node of the CN M 0 (M) mass of the considered chemical n (-) effective (i.e. secondary) rock formation porosity nc (-) number of clusters of particles having different travel times Q w (L 3 /t) pumping (or injection) well flow rate r i (L) pumping (or injection) influence radius r w (L) borehole radius Re (-) Reynolds number R h (L) hydraulic radius of the fracture flow Rρ (Ω⋅m) electrical resistivity of rock 0.1 0.9 , t t (t) time percentiles from pollutants breakthrough curves TD (-) time delay coefficient during pollutant mass-transfer U (L/t) velocity vector of the water flow x ij spatial oriented coordinate of fluid flow between nodes i and j of CNs X 1,1 (L 2 ) variance of one pollutant particle displacement including advection and pore-scale dispersion

Greek symbols
Δϕ (L) watertable drawdown (or mound) ϕ (L) water-pressure head Δm k (M) mass fraction of pollutant particles with same pathways  (t -1 ) rate coefficient of inactivation or biodegradation (or decay) of pollutant (or pathogen) λ d (t -1 ) scale invariant coefficient of pollutant mass transfer ν (L 2 /t) fluid kinematic viscosity   / (L⋅t) fluid viscosity/specific weight ratio  2 K (L 2 /t 2 ) variance of the hydraulic fracture conductivities  2 t (t 2 ) second moment of the expected concentrations in fractures; τ (-) tortuosity factor (<1) of the flow pathway in the fracture τ k (t) travel time of a generic particle cluster