Modeling radial groundwater flow in fractured media using fracture continuum approach

Two modeling approaches are commonly utilized for simulating flow in fractured formations: the discrete fracture network (DFN) approach and the stochastic continuum (SC) approach. Although the DFN approach is the most accurate, it has computational and memory constraints. The SC approach ensures fast processing but results in system over-homogenization. The fracture continuum (FC) approach arises as an integrated technique that incorporates the merits of both approaches. The main objective of this research is to develop a computationally efficient technique based on the FC approach to simulate the radial groundwater flow towards wells through two-dimensional fractured media under both steady and transient conditions. A stochastic generation of the DFN is performed in a Monte Carlo framework taking into account wells positioning. The DFN flow system is solved by applying the mass balance equation at fracture intersections. Fracture segments are mapped onto grids of 1 × 1 m and 5 × 5 m resolution as conductivity and specific storage cells. The grid flow problem is solved via MODFLOW. Flow and head discrepancies between the proposed technique and the DFN approach (reference solution) are assessed in steady and transient conditions. A grid-conductivity correction is needed to preserve the DFN flow in the presence of wells. A porosity estimation is proposed to identify the grid-pressure transient response. Promising flow and head results are observed for fine and coarse grid models. Some of the studied cases show large discrepancies in the maximum drawdown obtained on the coarse grid model. Accordingly, a new technique is proposed to handle such discrepancies and is found efficient in transient simulations (e.g., 11% and 26.12% discrepancies are minimized to − 0.93% and − 1.03% for two studied cases). The adopted mapping technique is found efficient when the interest is to estimate the average drawdown over an aquifer as correlation coefficients of 0.89 and 0.97 are found for the coarse and fine grid models, respectively when compared to the DFN model. However, the technique has limitations in estimating the drawdown at locations of wells.


Introduction
In a fracture-dominant medium, groundwater flow is mainly channeled through fractures. As the hydraulic behavior of a fractured rock is mainly based on the fracture system geometry, the in situ investigations of flow through fractured media showed that modeling the rock as a homogeneous continuum can create serious simplifications. Owing to the heterogeneities presented by fractures in a fractured rock mass, porous media-based modeling approaches have serious limitations. Consequently, a detailed description of the geometrical and hydraulic characteristics of the fracture system is needed to assess and understand the hydraulic behavior of fractured rock (Mo et al. 1998). These characteristics include length, orientation, density, connectivity, and permeability of fractures (Botros et al. 2008).
Responsible Editor: Broder J. Merkel Groundwater flow and contaminant transport in fractured systems have been active areas of research around the world for more than three decades now. Due to the nature of fractured systems, traditional porous media flow analysis methods are not directly applicable. Two different philosophies have been considered for simulating flow in fractured geologic media. These approaches are the discrete fracture network (DFN) approach and the stochastic continuum (SC) approach. The DFN modeling approach depends on the premise that groundwater flow occurs mainly within fractures. Darcy law is applied for each fracture segment and the mass balance must be achieved at each fracture intersection. The DFN approach was studied by Cacas et al. (1990), Sahimi (1995), Adler and Thovert (1999), Liu and Bodvarsson (2001), Park et al. (2001), Zhang and Sanderson (2002), de Dreuzy et al. (2002, 2004, Hyman et al. (2015), Maillot et al. (2016), Ahmed et al. (2019), Frampton et al. (2019), Jarrahi et al. 2019, Khafagy et al. (2020, Sweeney and Hyman (2020), and Khafagy et al. (2021). On the other hand, the SC modeling approach depends on the representative elementary volume (REV) concept which represents the rock mass as an equivalent porous medium. The hydraulic conductivities of the system are assigned to grid blocks based on random distribution. The SC approach was adopted in the studies of Neuman (1987Neuman ( , 1988, Neuman and Depner (1988), Follin and Thunvik (1994), Tsang et al. (1996), Niemi et al. (2000), Ohman and Niemi (2003), and Ando et al. (2003).
Each of the DFN and the SC approach has advantages and limitations. Although the DFN approach is more accurate in simulating the flow through fractured media, it is computationally very demanding. The characteristics of fractures must be incorporated into the DFN model as input parameters. The mass balance equation is applied at fracture intersections resulting in a huge number of linear equations. Thus, the DFN approach is not practical for simulating the flow in a large-scale problem. On the other hand, the SC approach is characterized by its processing speed, as it does not require detailed site-specific data about the fracture network (Ando et al. 2003). However, the SC approach over homogenizes the fractured rock. Consequently, this leads to large errors in characterizing flow-and transport-related quantities.
The fracture continuum (FC) approach that combines the advantages of the DFN and SC approaches and avoids their limitations has been developed in the early 2000s (Svensson 2001a;Svensson 2001b;Botros et al. 2008;Parashar andReeves (2009, 2011) and used recently by Parashar and Reeves (2011), Kalinina et al. (2014), Hadgu et al. (2015), and Ahmed et al. (2019). The FC technique depends on mapping the individual fractures of the network onto grid cells. Each grid cell is assigned an effective hydraulic conductivity in both directions and the system anisotropy is identified. The total flow is preserved for different cell sizes using a flow correction factor. The FC technique is an efficient, accurate, and fast processing tool. Nevertheless, it should be tested for other flow regimes such as radial flow (Botros et al. 2008).
The radial flow is one of the widely existing flow regimes in a subsurface flow system. Most studies of fractured aquifers had a great interest in the radial flow case, as it is economically important in water resources evaluation (Önder, 1998). The radial flow in fractured reservoirs was analyzed by Warren and Root (1963), Odeh (1965), and Kazemi (1969). Besides, the transient flow in fractured reservoirs was studied analytically. The radial flow was modeled stochastically through fractured porous formations by performing Monte Carlo simulations (Xu et al. 2013). Acuna and Yortsos (1995) presented the mass balance equation to solve the flow problem in discrete fracture networks in presence of wells and determine the groundwater heads at fracture junctions over the whole domain for both steady and transient conditions. In addition, a model describing the pressure transient in naturally fractured reservoirs, taking into consideration the effects of block size variations, was presented by Montazeri et al. (2011). Egya et al. (2018) studied the influence of variations in fractures' permeability and wells' locations in the discrete fracture network (DFN) on the pressure transient response. Furthermore, a pressure transient analysis was performed for complex discrete fractures utilizing semi-analytical models (Wang et al. 2018).
Addressing radial flow in a fractured system is crucial for many areas where fractured aquifers are abundant. Modeling such karst aquifer using traditional techniques has proven to be inaccurate and failed to reproduce aquifer responses to pumping tests.
Thus, the main objective of this study is to develop an efficient technique based on the FC approach to solve the radial flow problem in two-dimensional fractured media in presence of a single well or a system of wells. Also, addressing this issue for both steady and transient conditions and for different grid sizes is considered in this study so as to provide insight into the FC approach capabilities and how accurately or lack thereof, it can simulate the flow in these conditions. It is essential to emphasize that all the objectives will be studied in a stochastic Monte Carlo framework where different realizations are simulated using a Monte Carlo analysis. Different statistics of the fractured media will also be studied to assess the applicability of the developed model to a wide range of different fracture characteristics (fracture orientation, length, and permeability as well as density and connectivity of fractures).

Methodology
The adopted methodology to achieve the objectives of this research relies on the FC approach developed for natural groundwater flow conditions as in Svensson (2001a, b), Botros et al. (2008), Parashar andReeves (2011), andAhmed et al. (2019). This approach is modified here to account for stressed aquifers that are subject to pumping from a single well or multiple wells. The methodology implemented here can be summarized in several main steps. The first step is the random generation of the fracture network in a stochastic Monte Carlo framework based on fracture characterizations and locations of wells. From the created fracture network, the backbone fracture network (a connected network that does not contain dead-end fractures or closed loops) is obtained. In the second step, the flow system is solved on the DFN in both steady and transient conditions by applying the mass balance equation at fracture intersections. These DFN flow results are used as the reference solution to which the FC solution is to be compared. Fracture segments are mapped onto grids of 1 × 1 m and 5 × 5 m grids in the third step. In this step also, the discrete fractures are mapped onto cells to which hydraulic conductivity, as well as specific storage, are developed in a way that preserves flow characteristics. The fourth step involves the solution of the grid-based flow in steady and transient conditions via MODFLOW. In the fifth step, the grid flow and head results obtained on the mapped grid are compared to that of the DFN approach and the discrepancies are assessed in steady and transient conditions. In the last step, drawdowns are assessed over realizations and a new technique is proposed to minimize the errors of the maximum grid-based drawdown. Figure 1 presents a flowchart for the adapted methodology.

DFN generation in presence of wells
Modeling discrete fracture networks (DFNs) is mainly based on determining the site-specific data which is tough to obtain compared to porous media data. The fracture network geometrical data (e.g., fracture location, length, orientation, aperture, and spacing) are identified from the analysis of the outcrops, core samples, and aerial photographs. Also, geometrical data can be identified based on geophysical techniques such as ground-penetrating radar, boreholes imaging as well as seismic and electrical methods. The hydraulic properties of the DFN (e.g., transmissivity and specific storage) can be determined from hydraulic and tracer tests (Berkowitz, 2002). As the fracture properties vary spatially within the DFN and due to the uncertainty presented in the data, statistical distributions of these properties are developed. Many fracture networks are generated in a stochastic Monte Carlo framework based on the available data. In this study, the characteristics of fractures are assumed to be known and identified from field-specific data. Once they are known, the DFN generation starts.
Two intersecting fracture sets (perpendicular to each other) are generated within a two-dimensional synthetic domain. It is considered that uniform distribution is suitable for generating fracture starting locations (Berkowitz and Scher 1997). The probability density function (PDF) of the uniform distribution can be formulated as follows: where x is a random variant, a and b are the domain boundaries. While generating the x-coordinate of the fracture starting position in the synthetic domain, a and b are the minimum and maximum x-coordinates of the domain, respectively. The same is applied in the y-direction to generate the y-coordinate of the fracture starting position. As common in hydrogeological studies, fracture transmissivity is considered to follow a lognormal distribution (Botros et al. 2008;Ahmed et al. 2019). The PDF of the lognormal distribution can be formulated as follows: where T f is the fracture transmissivity (the random variable), ω and σ are the mean and standard deviation of the transmissivity logarithm, respectively. Each fracture set has a mean and a standard deviation of fracture log-transmissivities. Transmissivity does not change along the same fracture. In contrast, each individual fracture has distinctive transmissivity in the same set according to the lognormal distribution unless no variance is presented. This would result in spatial variations in the transmissivity of the fracture system. No flow interaction between rock matrix and fractures is considered. Fracture aperture can be calculated from the Cubic law which is a solution of the Navier-Stokes equation for flow between two parallel plates (Snow, 1969). The cubic law is defined as follows: where is the fluid density, g is the acceleration of gravity, μ is the fluid dynamic viscosity, and b is the fracture aperture. After sampling the fracture transmissivity from the lognormal distribution, it becomes easy to obtain the fracture aperture. In this study, the fracture aperture is considered to be incorporated into the fracture transmissivity distribution. In other words, the transmissivity distribution is used instead of the aperture distribution as a poor correlation is noticed between the mechanical aperture and the aperture presented in the cubic law (Bandis et al. 1985;Smith et al. 1987;Tsang, 1992). Fracture orientation is considered to follow a Fisher distribution (Nguyen et al. 2001). The PDF of the Fisher distribution (Fisher, 1953) can be formulated as: where ∅ is the fracture deviation from the set mean orientation and k is the Fisher constant. Fracture orientations ( become very close to the set mean orientation and tend to be normally distributed with variance equals (1/k ) in case of high k-value. Contrary, fracture orientations deviate further from the set mean orientation while decreasing the k-value till k=0. Fracture orientations tend to be uniformly distributed at k =0 (Best and Fisher, 1979). Every single fracture would have a ∅-value sampled from the Fisher distribution.
To get the fracture orientation, simply add the ∅-value, whether it was positive or negative, to the fracture set mean orientation. Fracture length is considered to follow either a power law or a lognormal distribution (Sahimi, 1995;Davy, 1997, 1998;Dreuzy, Davy and Bour, 2002;Botros et al. 2008). The PDF of the lognormal distribution is shown in Eq.
(2), but apply it for the fracture length ( l ) instead of T . Additionally, and σ (in this case) are the mean and standard deviation of the length logarithm, respectively. While the PDF of the power law distribution can be defined as follows: where a is the exponent of the power law, and l min is the minimum fracture length of the fracture set. Short fractures dominate the network connectivity when the power law exponent has a high value (a > 3.0) because multiple fractures would have lengths close to l min . In contrast, long fractures dominate the network connectivity when the power law exponent has a low value (a < 1.0). At 1.0 < a < 3.0, the network connectivity is dominated by both short and long fractures Davy, 1997, 1998). It is noticeable that the power law exponent ( a ) controls the network connectivity while sampling the fracture length from the power law distribution. Similarly, the variance of the length logarithm [ 2

log(l)
] controls the network connectivity while sampling the fracture length from the lognormal distribution. As 2 log(l) increases, long fractures control the network connectivity.
In this study, the locations of wells are assumed to be known. In addition, every single well is considered to be positioned at the intersection of two fractures. Before starting the network generation process, these two fractures are generated first. These two fractures are assumed to be long fractures of about ten times the mean length of the fracture set (deterministic lengths) to guarantee their interconnection to the backbone network. Furthermore, these two fractures have a common deterministic starting position which is chosen to be at the well location. The process of DFN generation stops when a pre-specified fracture density is achieved. The fracture density (m/m 2 ) can be calculated by dividing the summation of fracture lengths by the domain area. Ultimately, any parts of the fractures outside the domain will be trimmed.
To represent the actual connectivity of the fracture network, the fracture segments not participating in the network min l a for l ≥ l min flow (such as fracture dead ends and dead clusters that are not connected to the network) must be removed. To remove such segments, a percolation algorithm is applied. This should be done before mapping the discrete fractures onto grids to ensure eliminating the enhanced connectivity on the grid. Thus, eliminating the grid flow overestimation (Botros et al. 2008).

Steady-state flow solution using DFN approach
After applying the percolation theory, the backbone network is identified and the network actual connectivity is presented (Robinson, 1983;Englman et al. 1983;Hestir and Long, 1990;Berkowitz and Balberg, 1993;Berkowitz, 1995;Botros et al. 2008;Ahmed et al. 2019). The steady-state flow problem is solved over the backbone network within a two-dimensional synthetic domain where the right and left boundaries are assumed to be specified head boundaries while the upper and lower boundaries are considered as noflow boundaries. However, the flow will not be in one main direction (as expected in porous media) on the microscopic scale due to the anisotropy and heterogeneity of fractures. The mass balance equation in steady conditions (Acuna and Yortsos, 1995) is applied at the fracture intersections of the backbone network.
where j indicates all joints connected to jointi , i is the fluid density at jointi , g ij is defined as the segment conductance between joints j and i and is calculated as ( g ij = A K ) , A denotes the cross-sectional area of the fracture and it equals to the fracture aperture (b) as the domain-third dimension is taken as unity, is the specific weight of fluid, K is the fracture hydraulic conductivity, (P j − P i ) denotes the pressure drop between joints j and i which is estimated as H j − H i , (H j − H i ) indicates the head difference between joints j andi , l ij is the fracture segment length between joints i and j , Q w is the well pumping rate, is defined as the Kronecker delta, and m is the joint number of the well. im = 1 .0 fori = m , otherwise im = 0.0. AsT = Kb , Eq. (6) can be simplified as follows: A large number of linear equations arises. Once the flow problem is solved and the head values at the fracture intersections are obtained, the flow of each fracture segment ( Q ij ) can be computed using Darcy's law as follows: The flow of the backbone network is solved over multiple realizations. The global flow from the left to the right boundary is calculated and is referred to as Q DFN . DFN flow results are used as the reference solution to which the FC solution is to be compared.

Transient flow solution using DFN approach
The transient groundwater simulations are required to predict the variations in water levels with time. The transient response of fracture systems can be described using the DFN technique. The transient flow problem is solved by applying the mass balance equation in the transient conditions at the fracture intersections of the backbone network (Acuna and Yortsos, 1995).
where c f indicates the fluid compressibility which is the reciprocal of the fluid bulk modulus, ΔP i is the difference between the pressure at the current time and that at the previous time at node i , Δt is the time step, and V i is the share of fluid volume in the segments to node i . The shared fluid volume ( Vi) can be computed as follows: The bulk modulus of water equals 2100 MPa (Terzaghi and Peck. 1967). Based on Eq. (7), Eq. (9) can be simplified as follows: The DFN model is firstly solved under steady-state flow conditions, then the head outputs are used as initial heads for the transient DFN model.

Mapping fractures onto grid
Fracture mapping techniques have been developed owing to the complexity and computational constraints of the DFN technique. This study depends mainly on Botros et al.'s (2008) proposed approach which differs from Svensson's (2001a) approach. Svensson's (2001a) approach has limitations related to the ratio of fracture aperture to the cell size these fractures are mapped on. Furthermore, it leads to different values of hydraulic conductivity at head cell faces intercepted by the same fracture which is not reasonable. Botros et al.'s (2008) approach is more comprehensive. It overcomes the limitations of Svensson's (2001a) approach. In addition, it takes into consideration the correction factors presented in the works of Langevin (2003) and McKenna and Reeves (2006). To map a discrete horizontal or vertical fracture on a finite-difference grid of cell size ( Δ) , the hydraulic conductivity (in the direction of fracture) of the cells intercepted by this fracture can be computed as follows: Mapping oriented fractures on a grid leads to a "stairstep" pattern. This pattern makes the flow path on grid longer than on DFN. Consequently, the head gradient is decreased. Thus, a correction factor is introduced in Eq. (13) to estimate the fracture contribution to the hydraulic conductivity of the cell intercepted by the inclined fracture. The correction factor compensates for the decrease in head gradient by an increase in the cell hydraulic conductivity so that the same amount of flow can be preserved (Botros et al. 2008;Parashar and Reeves, 2011). The fracture contribution to the grid cell conductivity (K b ) can be calculated as follows: where is the angle of fracture with the x-direction and C( ) is the correction factor used to estimate the cell hydraulic conductivity in case of mapping an oriented fracture. According to Botros et al.'s (2008) study, fractures are classified into two types (type I and type II fractures) depending on their alignment with the finite-difference grid cells. Type I fractures intercept two cell opposite faces. The contribution of such fractures to the hydraulic conductivities of the cell ( K b ) is in a single direction which is orthogonal to the intercepted cell faces. The other-direction hydraulic conductivity is taken as the matrix conductivity. On the other hand, type II fractures intercept two adjacent faces of the cell. In case of a single type II fracture crossing a cell, the contribution of such fracture to the cell hydraulic conductivities ( K b ) will be in both directions. Nevertheless, in case of multiple type II fractures are existing within a particular cell, the contribution of these fractures to the cell conductivities is divided between x-direction and y-direction equally.
To obtain the cell hydraulic conductivity in a specific direction, the cell hydraulic conductivities in that direction due to the interception of many fractures to the cell are accumulated. Therefore, cell hydraulic conductivity in a specific direction. Eventually, each cell would have a distinctive anisotropy ratio. It is necessary to emphasize that the pre-mentioned mapping rules are applied for mapping the discrete fracture network on a finite-difference grid and obtaining the domain conductivities.
The flow path of the fractures on a grid becomes representative of the actual flow path in case of mapping fracture zones of high density. The fracture cell is surrounded by cells of high conductivities in both directions. The original orientations of fractures in such zones are hard to investigate. Therefore, a loss in directionality arises. To identify the cell as a part of a dense fracture zone, two conditions have to be achieved. First, the cell ought to be part of a squared continuum area. The four cells of the squared area should be type II fracture cells where the hydraulic conductivities are assigned in both directions. Thus, Botros et al. (2008) introduced a modified correction factor C( ) as follows: Matrix transmissivity is assumed at least six orders of magnitude lower than the fracture set transmissivity to ensure a great contrast when compared to the fracture transmissivity. Thus, the flow is channeled through fractures. The global flow over the grid Q GRID is determined using MOD-FLOW, then compared to Q DFN . The flow overestimation correction factor is calculated as the average of Q GRID ∕Q DFN over realizations. The hydraulic conductivities of the cells are divided by this correction factor to preserve the system flow (Botros et al. 2008).

Steady and transient flow solutions over the grid
Once the mapping technique is applied and the DFN is converted to hydraulic conductivities on finite-difference grid cells after applying the flow overestimation correction factor, the partial differential groundwater flow equation can be solved over the grid in whether steady or transient conditions. The boundary conditions, locations of wells, pumping rates of wells, and conductivity distribution are assigned to the finite-difference grid. The partial differential equation of two-dimensional transient groundwater flow in a heterogeneous and anisotropic confined aquifer in presence of pumping wells can be described as shown in Eq. (16). This partial differential equation is solved numerically using MODFLOW (Harbaugh et al. 2000).
where h is the groundwater head, Q represents the pumping rate of the well (sink rate) at defined coordinates and specific time, and S s denotes the specific storage of the aquifer which can be calculated as follows: where is the aquifer compressibility, n e is the effective porosity, and C f is the fluid compressibility as stated earlier.
The aquifer compressibility as well as the specific yield of the aquifer are assumed to be negligible in this study. The steady and transient flow problems are solved over fine and coarse grid resolutions. The grid-based model is solved in steady conditions after applying the flow overestimation correction factor, then the output heads are assigned as initial heads to the transient grid-based model. In steady-state flow conditions, Eq. (16) can be written as follows: MODFLOW code is primarily written in Fortran programming language. Therefore, using Fortran to implement the proposed approach is efficient in achieving direct linkage to the MODFLOW code. For each realization, the developed code converts the geometrical and hydraulic parameters of the discrete fractures to grid distributions of hydraulic conductivity, anisotropy, and specific storage. Besides, these distributions, for different cell sizes, are generated in several output files which are later used as input files for MODFLOW code to solve the grid-based flow problem for the different grid resolutions in a sequential manner. Ultimately, the whole Fortran code (i.e., integration between the developed code and MOD-FLOW code) automatically repeats the same process for a predefined number of realizations.

Porosity estimation over grid
One of the most important hydraulic parameters in case of studying the transient flow response is the system porosity. The changes in the specific storage of the aquifer are mainly dependent on the aquifer porosity. Specific storage is assigned to fracture cells only as the matrix is considered impervious. Consequently, the specific storage is assigned as zero to the matrix cells. The volume of water stored in the discrete fractures is usually less than that stored in the fracture cells these fractures are mapped on. Thus, the total volume of water stored in the discrete fractures is distributed over the grid fracture cells. Ultimately, a porosity estimation over the fracture cells is proposed as follows: where ∑ V i represents the total volume of water stored in the discrete fractures and N FC indicates the total number of fracture cells. Δx , Δy , and Δz are the grid discretization in x-direction, y-direction, and z-direction, respectively. As stated earlier, the aquifer depth is considered as unity (i.e., Δz = 1.0).

Flow and head comparisons
After the DFN and the grid-based models are solved, the grid global flow of each realization can be computed and compared to that of the DFN model to show the errors. Furthermore, the average and standard deviation of the gridbased model flow over a large number of Monte Carlo realizations are calculated after applying the flow overestimation correction factor resulting from the solutions of a limited number of realizations, then compared to that of the DFN model. Also, a Fortran code is written to compare the head results of the DFN model with that of the grid-based model for the same realization.

Assessment of drawdown errors
Head and drawdown distributions are expected to change between realizations. Consequently, errors of head and drawdown must be assessed over the generated realizations. Errors are assessed qualitatively by comparing the average, standard deviation, and maximum drawdown of the gridbased model to those of the DFN model. In addition, histograms for the errors in the average drawdown over realizations are introduced as quantitative measures. Also, the maximum drawdowns of the transient grid-based model with time are compared to the ones of the transient DFN model of the same realization. Moreover, the same comparison is done for the averages and standard deviations of the system drawdowns with time. The percentage error in the average drawdown for one realization is computed as follows: where AvgDD DFN and AvgDD Grid are the average drawdown of the DFN model and the grid-based model, respectively.

Proposed technique to minimize the errors in the transient grid maximum drawdown
In the DFN model, the well extracts water from fractures that look like narrow pipes while in the grid-based model, the well extracts water from very wide grid cells (enhanced connectivity) compared to the aperture size. So, the impact of well extraction is expected to be higher on DFN than on grid. In other words, the drawdown of the DFN model at the well location would be larger than that of the grid-based model. Accordingly, a technique is proposed to minimize the error in maximum drawdown on grids. It depends on the numerical approximation of the two-dimensional flow equation for heterogeneous, anisotropic confined aquifer (Konikow, L.F., 2008) which is defined as follows: where T xx[i− 1 2 ,j] represents the average transmissivity between nodes i, j and i − 1, j , T xx[i+ 1 2 ,j] represents the average transmissivity between nodes i, j and i + 1, j . Similarly, for T yy[i,j− 1 2 ] and T yy[i,j+ 1 2 ] but in Y-direction. h i,j,n is the head value at node i, j at time stepn , h i−1,j,n is the head value at previous node in X-direction, h i+1,j,n is the head value at next node in X-direction, h i,j−1,n is the head value at previous node in Y-direction, h i,j+1,n is the head value at next node in Y-direction, and h i,j,n−1 is the head value at the same node but at the previous time step.
The proposed technique depends on obtaining the DFN maximum drawdown and DFN heads at the position of maximum drawdown in steady conditions before and after pumping is implemented. After the grid-based model is solved in steady conditions, the grid maximum drawdown is obtained and error is assessed. A target head ( h target )-in case of wells are on-is obtained by subtracting the DFN maximum drawdown from the grid head at the same location in case of no extractions. The ratio between h i−1,j,n − h i,j,n and h i−1,j,n − h target is calculated and the result is multiplied by T xx[i− 1 2 ,j] . Also, the ratio between h i+1,j,n − h i,j,n and h i+1,j,n − h target is multiplied by T xx[i− 1 2 ,j] . Ultimately, the hydraulic conductivities of the well cell and its two adjacent cells in the main flow direction (X-direction in this study) are estimated. The grid-based model is run in transient conditions and the grid maximum drawdown is finally obtained at different time steps. (21)

Results and discussion
Four different cases are studied in steady and transient conditions as shown in Table 1. The DFN parameters introduced in Table 1 are pre-defined in "DFN generation in presence of wells." Case 1 is studied in presence of a single pumping well at the middle of the domain (scenario 1) as well as in presence of five pumping wells (scenario 2). All other cases are studied in presence of five pumping wells only. Locations of the pumping wells are presented in Table 2. It is noticeable that the wells are allocated symmetrically to examine the influence of the heterogeneity and anisotropy of fractures on the drawdown distribution.

Mapping results
The backbone fracture network is mapped to a two-dimensional grid with cell sizes 1 × 1 m and 5 × 5 m for all cases introduced in Table 1. To simulate transmissivity contrast between fracture cells (with high flow velocity) and matrix cells (with no or very low flow velocity), the matrix cells' transmissivity is assumed to be six orders of magnitude lower than the mean transmissivity of fractures (Botros et al. 2008). The distributions of grid conductivities for one realization of case 3 on cell sizes 1 × 1 m and 5 × 5 m are shown in Fig. 2 as a sample. It could be noticed that mapping discrete fractures on coarser grid resolutions would add more homogenization to the system. In addition, the system conductivity on the coarse grid was found smaller than that on the fine grid so that the system flow could be preserved. Spatial variations in the grid conductivities and conductivity  contrast between the matrix and fracture cells could also be observed.
After obtaining the porosity distribution as presented in "Porosity estimation over grid," the specific storage distribution was directly obtained. Figure 3 shows the distribution of specific storage over a grid of cell sizes 1 × 1 m and 5 × 5 m for a single realization of case 1. It could be seen that the specific storages are assigned to fracture cells only while its values are assigned as zeros to matrix cells. Moreover, the specific storage values over the grid of cell size 5 × 5 m are smaller than those over the grid of cell size 1 × 1 m because the volume of water is distributed over a larger area in the case of 5 × 5 m grid resolution resulting in smaller grid porosity values.

Verification results in steady conditions
The flow problem of 100 realizations was solved over the grid in the steady-state flow condition for each studied case using MODFLOW. The grid-based model solutions were verified by comparison with the results of the DFN model. The pumping rate of the single well (W1) was assumed to be 250 m 3 /h for scenario 1 while in the other cases, where five wells are existing, the five wells were assumed to have equal pumping rates of 50 m 3 /h each. Once the head distribution was determined over the DFN as well as over the grid for each realization, the global flows of the DFN model ( Q GRID ) and the grid-based model ( Q DFN ) could be calculated by applying Darcy's law. The flow overestimation factor (the conductivity correction factor) was estimated for each grid resolution of the four cases as the average of Q GRID ∕Q DFN for the 100 realizations as shown in Table 3. Larger values for flow overestimation factor were observed for the coarser grid resolution due to the enhanced connectivity of the mapped fractures. Figure 4 shows comparisons between total grid flows ( Q GRID ) at the pre-mentioned two grid resolutions and the total network flows ( Q DFN ) for 100 realizations for case 1 (scenario 1), case 1 (scenario 2), and case 3, respectively. It could be seen that, for all cases, a good agreement was observed between Q DFN and Q GRID for the fine grid resolution while a flow overestimation was observed for the coarse grid resolution because of the larger enhanced connectivity that is present in coarser fracture cells. The flow results had slight changes for the two scenarios of case 1. This ensures that the technique used for generating the fracture network (based on locations of wells) does not add large discrepancies to the flow results. Consequently, the proposed approach is robust while solving the flow in presence of either a single well or multiple wells.   In addition, slight variations in the values of the flow overestimation factor between realizations were observed. Therefore, according to Botros et al. 2008, a larger number of realizations (200 realizations in this study) was generated and solved over grids with adjusted hydraulic conductivity based on the flow overestimation factor of the 100 realizations. Figure 5 shows comparisons between the grid total flow and the network total flow for both the fine and coarse grid resolutions after adjusting the grid conductivity for the same cases previously presented in Fig. 4. It could also be noticed from Fig. 5 that the flow overestimation was minimized, and a good agreement was observed between the grid flow results and the DFN flow results.
The DFN and the grid-based models were solved. The head results of the grid-based model were compared to those of the DFN model for the same realization of the tested case. Head values of the grid-based model were assigned to the joints of the DFN, then these grid heads at the positions of the backbone joints were compared to the heads of the DFN model. The comparisons between the heads of the grid-based model (for 1 × 1 m and 5 × 5 m grid resolutions) and those of the DFN model for a single realization (the same realization for both resolutions) of case 3 are shown in Fig. 6. A strong correlation was found between the DFN heads and the grid heads. Mapping fractures onto coarse grid made the grid head points more dispersed around the  45° line which had been used to test the equivalence between grid and DFN heads.
The drawdown distributions for a single realization (Reaz1) of the two scenarios of case 1 as well as for another realization (Reaz2) of scenario 2 of case 1 are shown in Fig. 7. Unlike flow simulations in porous media, Fig. 7 shows that the drawdown is not distributed symmetrically over the domain. This is mainly due to the heterogeneity and anisotropy that fractures bring to the system. Although no significant change was observed between the flow results for the two scenarios of case 1, a significant change in the drawdown distributions of the two scenarios was found. The maximum drawdown of the first scenario was 1.44 m while it was 0.57 m for the second scenario. In addition, the maximum drawdown of the second scenario (in case of five wells are existing) is not definitely at the middle of the domain unlike the case of simulating flow through homogeneous and isotropic porous media. This may happen in some realizations (as shown in Fig. 7) due to system heterogeneity and anisotropy and may also happen due to transmissivity variability. Although the maximum drawdown of one realization (Reaz1) of scenario 2 of case 1 is at the middle of the domain, it is at (250 m, 750 m) for another realization (Reaz2) as shown in Fig. 7.

Verification results in transient conditions
The DFN model and the grid-based model were solved in presence of pumping wells for single realization in steadystate flow conditions. Then, the conductivity correction factor for the realization was estimated. It was utilized to adjust the grid conductivities of the transient flow model. The output heads of a steady-state model in presence of no wells were assigned as initial heads to the transient flow model. Moreover, the domain-specific storage was computed and assigned to the model after estimating the grid porosities. The transient flow and drawdown results of the grid-based model were compared to those of the DFN model. All cases introduced in Table 1 were studied under transient conditions. The mean transmissivity of the fracture sets was assumed to be 10 −9 m 2 /s for all cases presented in Table 1. It is important to clarify that, in this synthetic model, the value of the mean transmissivity was chosen to be three orders of magnitude less than that presented by Botros et al. (2008) because the hydraulic gradient was very high (i.e., 0.03). The combination of this hydraulic gradient and the mean fracture transmissivity presented by Botros et al. (2008) caused the system to reach steady-state conditions in a very short time. This attributes to the resulting relatively high inflow values and the rapid response of fractures. Consequently, the pre-mentioned mean transmissivity (i.e., 10 −9 m 2 /s) is used to express the system transient response of the DFN and the grid-based models. The pumping rate of the single pumping well (W1) was assumed to be 10 −9 m 3 /s in case 1 while in the other cases, where five wells are existing, the five wells were assumed to have equal pumping rates of 2 × 10 −10 m 3 /s each. For case 1, the simulation time was considered as 5 days and the time increment was chosen as 2 h while the simulation time was considered as 2 days and the time increment was chosen as 1 h for other cases. Such small values of pumping rates were chosen to be consistent with the small values of the mean fracture transmissivity in order to study the transient response of the system because higher values of pumping rates lead to a very fast change from transient to steady-state condition.
The conductivity correction factor was applied for the steady-state model in case of no wells. The output heads of the steady-state model were assigned as initial heads to the transient flow model. Figure 8 shows a comparison between the initial heads (at t = zero) of the DFN transient model of case 1 (scenario 1) and those of the grid-based transient model for grid resolution 1 × 1 m. A good agreement was observed between the initial heads of both models. Moreover, the conductivity correction factor was applied to the  (19,17). The results of the transient flows (inflow and outflow of the system) of the grid-based model (for 1 × 1 m and 5 × 5 m grid resolutions) and the DFN model for scenario 1 of case 1 were compared. Such comparisons were shown in Fig. 9. It could be noted that the system inflow increased with time while the system outflow decreased. This resulted from the increased drawdown with time while the extraction was constant. The increased drawdown caused an increase in the hydraulic gradient between the left boundary and the well. Also, it caused a decrease in the hydraulic gradient between the right boundary and the well. A very good agreement was found between the transient flow results of the DFN model and the grid-based model.
The drawdown distribution of the grid-based transient model was compared to that of the DFN transient model for case 1 (scenario 1) for 1 × 1 m grid resolution. Figure 10 shows such comparisons at times equal to 0.167, 0.75, and 5 (end of simulation time) days, respectively. Strong correlations were found between the drawdown distributions of both types of models at the times chosen for comparison. To ensure the agreement between the grid-based transient model and the DFN transient model at all times, some statistical analyses were performed to show the drawdown distributions of both models with time. These analyses were the maximum drawdown with time, the average drawdown with time, and the standard deviation in the domain drawdown with time. Such analyses are presented for 1 × 1 m and 5 × 5 m grid resolutions for case 1 (scenario 1) in Fig. 11, respectively. Seemingly, errors in the grid drawdown results for case 1 (scenario 1) were not significant for both grid resolutions. This mainly happened because there was no variability in the fracture transmissivity. This will be shown in Table 4 for the maximum grid drawdown results.
It is essential to emphasize that the average grid drawdown, as well as the standard deviation of the grid drawdown, are computed after excluding the drawdown values of the matrix cells. Good agreements were observed between the grid-based model and the DFN model regarding the maximum, average, and standard deviation of drawdown with time. Furthermore, the drawdown increased exponentially with time till reached a steady-state condition. Larger errors were observed between the maximum drawdown with time for the 5 × 5 m grid resolution model and that of the DFN model compared to the errors noticed for the 1 × 1 m grid resolution model. The main reason is that more homogeneity was added to the system in case of the 5 × 5 m grid resolution model compared to the 1 × 1 m grid resolution model. In addition, in case of the 5 × 5 m grid resolution model, the head was distributed over larger cell sizes. Table 4 present the maximum grid drawdown for 1 × 1 m and 5 × 5 m grid resolutions, the maximum DFN drawdown, and the absolute percentage errors in the maximum grid drawdown for the test cases. The results of the maximum drawdown showed good agreement for the fine grid resolution. On the other hand, large errors arose in the maximum grid drawdown in the case of coarse grid resolutions especially in cases 3 and 4 where the transmissivity variability was presented. In general, coarser grid models showed greater discrepancies in the maximum

Assessment results of drawdown errors
Drawdown distributions were found changeable between realizations. As such, drawdown errors were computed over 100 realizations and assessed qualitatively and quantitively. The average and standard deviation of the grid drawdown for 100 realizations of case 3 for 1 × 1 m and 5 × 5 m grid resolutions were compared to those of the DFN drawdown as shown in Fig. 12, respectively. The correlation coefficient of the average drawdown of 1 × 1 m and 5 × 5 m grid resolution models were 0.97 and 0.89 (strong correlations), respectively. The correlation coefficient of the standard deviation of the drawdown of 1 × 1 m and 5 × 5 m grid resolution models were 0.92 and 0.86 (strong correlations), respectively. Lower correlations (i.e., larger errors in the average grid drawdown) present in the coarser grid models were ascribed to their accompanying enhanced fracture connectivity in addition to local errors at locations of wells (maximum drawdown) as described in Table 4. The maximum drawdowns and the heads at maximum drawdown locations of the grid-based model (for 1 × 1 m and 5 × 5 m grid resolutions) were compared to those of the DFN model in case of wells are off as well as in case of wells are on for 100 realizations as shown in Fig. 13. Strong correlations were observed for the fine grid model (1 × 1 m).
In contrast, relatively moderate correlations (moderate to strong correlations) were found for the coarse grid model (5 × 5 m) due to the enhanced connectivity of the mapped fractures on the coarse grid. Ultimately, a great similarity was found between the errors of the heads at maximum drawdown in case of wells are off and those in case of wells are on which means that adding wells to the system does not add errors to the head results at the maximum drawdown location.
Histograms for the errors in the average drawdown over the 100 realizations for 1 × 1 m and 5 × 5 m grid resolutions are shown in Fig. 14. About 98% of the realizations gave errors in the average drawdowns between − 10% and + 10% for 1 × 1 m grid resolution, while 83% of the realizations gave errors within the same error range for 5 × 5 m grid resolution. The average percentage error over the 100 realizations was 2.59% for 1 × 1 m grid resolution, while it was 5.86% for 5 × 5 m grid resolution. Obviously, the fine grid models gave better results than the coarse grid models due to the enhanced connectivity accompanying the coarse grid resolutions as the hydraulic properties of fractures are mapped over larger areas. Comparison between the grid drawdown (for 5 × 5 m grid resolution) and DFN drawdown for a single realization of (a) case 3 and (b) case 4 before and after correction at end of simulation time

Results of the proposed technique to minimize the errors in the transient grid maximum drawdown
Errors in the maximum grid drawdown were observed. Thus, a technique was proposed to minimize such errors. As presented in Table 4, the errors in the maximum drawdown of a single realization of case 3 and case 4 reached 11% and 26.12% at end of simulation time (when the flow became steady) for 5 × 5 m grid resolution. The proposed technique for minimizing the error in maximum drawdown was implemented for the pre-mentioned realizations of cases 3 and 4. Figure 15 shows a comparison between the grid drawdown (for 5 × 5 m grid resolution) and the DFN drawdown before and after applying the proposed technique for a single realization of case 3 and case 4, respectively. The errors in the maximum drawdown of the realizations of cases 3 and 4 were minimized to − 0.93% and − 1.03% at end of simulation time. The domain drawdown (except at the well location) was not affected significantly.
It essential to emphasize that the proposed technique is not fully developed and needs further development because it depends on having the maximum DFN drawdown. However, it was useful in transient simulations as the DFN model could be run in steady conditions and the maximum drawdown could be obtained. The corrected conductivities of the cells could be assigned in the transient grid-based model. The total runtime of the steady-state DFN model and the transient grid-based model would be relatively much less than the runtime of the transient DFN model. Figure 16 shows a comparison between the maximum drawdown with time of the 5 × 5 m grid-based model and that of the DFN model before and after applying the proposed technique for one realization of case 3 and case 4, respectively. A good agreement was found between the maximum grid drawdown and the maximum DFN drawdown after correction. Furthermore, the proposed technique may be a step to start a transport grid-based model.

Summary and conclusions
This study intends to develop a computationally efficient technique based on the fracture continuum (FC) approach to model the radial groundwater flow through two-dimensional fractured media in steady and transient conditions. A generation of the DFN is proceeded in a Monte Carlo framework taking into consideration the locations of wells. The DFN flow problem is solved in steady and transient conditions. The DFN is mapped as hydraulic conductivity and specific storage on cells of sizes 1 × 1 m and 5 × 5 m. The grid flow problem is solved via MODFLOW in steady and transient conditions, then verified against the DFN solution, and the upscaling technique's discrepancies are assessed.
Promising flow and drawdown results are observed over 1 × 1 m grid resolution models (fine models). On the other hand, large flow and maximum grid drawdown errors are found over 5 × 5 m grid resolution models (coarse models). A grid-conductivity correction is needed to preserve the DFN flow in presence of wells. Although the wells are symmetrically assigned, the drawdown resulting from the proposed technique is asymmetrically distributed, unlike simulating fracture systems as homogeneous media (REV concept). Moreover, a porosity estimation is proposed to present the pressure transient response over the grid system. Drawdowns are assessed over realizations. In some cases, large errors are found in the drawdown results of the coarse grid models at locations of wells. Accordingly, a three-cell conductivity correction technique is proposed and promising results regarding the maximum grid drawdown are noticed (e.g., 11% and 26.12% errors are minimized to − 0.93% and − 1.03% for two studied cases). This technique needs further development, but it is found efficient in case of transient simulations as it saves significant runtime. Ultimately, the adopted mapping technique is found very efficient when the interest is to estimate the average drawdown over an aquifer as correlation coefficients of 0.89 and 0.97 are found, for one of the studied cases, for the coarse and fine grid models, respectively, when compared to the DFN model. Moreover, 98% and 83% of the realizations give errors in the average drawdowns between − 10% and + 10% for the fine and coarse grid resolutions, respectively, and the corresponding average percentage errors over the 100 realizations are 2.59% and 5.86%. Contrarily, the approach has limitations in estimating the maximum drawdown or the drawdown at locations of wells.
For related future perspectives, it is recommended to investigate the mapping technique in case of different boundary conditions, study the grid refinement around wells (i.e., mapping discrete fractures on grid of variable cell sizes), consider the matrix-fracture flow interaction, and extend the proposed technique to model the transport in fractured media in presence of wells.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.