A dual-continuum model (TOUGH2) for characterizing flow and discharge in a mechanically disrupted sandstone overburden

Underground hard coal mining usually disrupts the mechanical equilibrium of rock sequences, creating fractures within minor permeable rocks. The present study employs a dual-continuum model to assess how both fractured and porous sandstone media influence the percolation process in postmining setups. To test the approach, the software TOUGH2 was employed to simulate laminar fluid flow in the unsaturated zone of the Ibbenbüren Westfield mining area. Compared to other coal mining districts in Germany, this area is delineated by the topography and local geology, leading to a well-defined hydrogeological framework. Results reveal good agreement between the calculated and measured mine water discharge for the years 2008 and 2017. The constructed model was capable of reproducing the bimodal flow behavior of the adit by coupling a permeable fractured continuum with a low-conductivity rock matrix. While flow from the fractured continuum results in intense discharge events during winter months, the rock matrix determines a smooth discharge limb in summer. The study also evaluates the influence of individual and combined model parameters affecting the simulated curve. A detailed sensitivity analysis displayed the absolute and relative permeability function parameters of both continua among the most susceptible variables. However, a strong a priori knowledge of the value ranges for the matrix continuum helps to reduce the model ambiguity. This allowed for calibration of some of the fractured medium parameters for which sparse or variable data were available. However, the inclusion of the transport component and acquisition of more site-specific data is recommended to reduce their uncertainty.


Introduction
Underground hard coal mining operations irreversibly disrupt the preexisting mechanical equilibrium of the geological media (Kim et al. 1997;Newman et al. 2017). The employment of high-recovery methods redistributes, concentrates and reorientates the stress state of the bedrock above the mined seams (Palchik 2003;Meng et al. 2016;David et al. 2017). Conceptually, coal mine overburden can be subdivided into caved, fractured and deformation zones ( Fig. 1), according to their response to the disruptions (Palchik 2003;Zhang et al. 2018a;Bai and Tu 2019;Liu et al. 2019). The caved zone encompasses the highly fragmented area formed after the overburden falls to the mine level. Above it, horizontal and subvertical fractures break the bedrock, developing a fractured zone, whereas towards the top of the rock sequence, layers of higher plasticity create a deformation zone, dampening the movement of the underlying sequence. These overburden deformations also modify the hydrogeological setting of mined areas. Induced fractures, for instance, alter the flow paths of percolating groundwater while increasing the hydraulic conductivity of the rocks (Qu et al. 2015;Zhang et al. 2018b;Liu et al. 2019;Qiao et al. 2019). In addition, subsidence related to the deformation zone may change surface hydrological networks, increasing or decreasing the water balance of the mined area and changing the water level in the aquifers (Kelleher et al. 1991;Blodgett and Kuipers 2002). Accordingly, fluid flow in postmining overburdens can be described by four processes: (1) meteoric water infiltration into the first meters of soil and disaggregated rock, (2) percolation through an undisturbed zone, (3) groundwater flow in subvertical fracture networks, and (4) flow and groundwater storage in the porous matrix (Fig. 1). However, as the latter two processes take place in the same physical space, they are rarely modeled simultaneously. The simulation of a discrete setup that considers fractures and porous matrix requires extensive characterization of both media (Ghasemizadeh et al. 2012;Hao et al. 2013). This results in models with poor real-field practicability and lots of uncertainties due to the high demand for restricted-access data, especially those related to fracture size, spacing, orientation and transmissivity distributions (Kovács and Sauter 2007;Koohbor et al. 2019). Addressing the influence of such uncertainties and, therefore, their effect on the models' predictive capability is a complex task that requires the use of probability distribution methods such as Monte Carlo or polynomial chaos (Rouleau and Gale 1987;Younes et al. 2020;Drouin et al. 2021;Guo et al. 2021). In the presence of a large number of variables, the model requires a significant amount of simulations to cover the range of parameters, being computationally robust and expensive. Furthermore, if poor parameters and constraints are introduced into the model, poor results will be obtained (Fielding et al. 2011). All this has led most postmining models to explicitly ignore fracture structures by employing the equivalent porous medium (EPM) approach or even both media with lumped parameter models (e.g., see Adams and Younger 2001;Banks 2001;Vincenzi et al. 2010). Notwithstanding that both approaches have been successfully applied in case studies, they exhibit notable drawbacks when describing the fast component of the fractured system (Kim et al. 1997;Rapantova et al. 2007).
The present study proposes the dual-continuum (DC) approach as a potential tool to include fracture and matrix flow in minewater models. The DC approach was first proposed by Barenblatt et al. (1960) to describe fluid seepage  David et al. 2017 andBai andTu 2019) in fissured rocks without treating the subsurface as a single porous medium. The approach overcomes the difficulties of obtaining detailed information for constructing discrete models by handling fractures as a network of averaged characteristics. Thus, fractures and porous media are considered as two overlapping and interacting continua, with different flow, transport and storage parameters. The interaction between the two is achieved through a mass transfer function determined by the size and shape of the blocks, as well as their local difference in pressure, temperature and chemical potentials (Pruess and Narasimhan 1985;Beyer and Mohrlok 2006;Aguilar-López et al. 2020). Over time, the DC approach has been adapted for numerous subsurface processes, including oil recovery, geothermal energy, nuclear waste repositories, and CO 2 sequestration (e.g., see Kasiri and Bashiri 2011;Wu et al. 2011;Azom and Javadpour 2012;Hao et al. 2013). However, it has been particularly suited for the simulation of diffusive flow and rapid drainage processes in two-dimensional (2D) karst aquifers (Sauter 1992;Kovács and Sauter 2007;Kordilla et al. 2012;Dal Soglio et al. 2020).
Here, the software package TOUGH2 is employed to test the applicability of the DC approach in mining-derived flow regimes characterized by induced fractures within porous units. For this purpose, the well-defined hydrogeological system of the Ibbenbüren Westfield in north-west Germany is selected. Although the area was mined for decades within sharp topographic and geologic boundaries, only a few hydrogeological models have been published. By the time the Westfield mine was active, reports were limited to pointing out fluid flow through a fractured sequence due to low measured permeabilities, below 1 mD (e.g., see Lotze et al. 1962;Bässler 1970). The first reference to a numerical model appeared for the mine closure plan in 1979, where water rebound was calculated using boxes to accommodate the water within empty spaces left by mining. However, this model was questioned some time after, as the rebound occurred 2 years later than predicted (Klinger et al. 2019). Some of the reasons given for the delay were the disregard of water-bearing fractures along with the imbibition process of unsaturated permeable units. Most recently, Rudakov et al. (2014) presented a lumped parameter model to simulate the temporal distribution of the mine water discharge as a factor of the precipitation. This model, however, shows deviations of up to 40% in summer months for heavy precipitation events. Lastly, Coldewey et al. (2018) developed an equivalent continuous model with a fixed recharge boundary to replicate the measured phreatic points around Westfield. This fit was achieved by setting rock permeabilities around 1 × 10 -12 m 2 , which are values three orders of magnitude above the maximums reported in the literature. With this hindsight in mind, the present DC model configures a breakthrough in the application of flow modeling by considering and quantifying the specific contribution of fractured and matrix continua to the final transient mine discharge signal of the entire coal basin. Likewise, the study further evaluates the effect of individual and combined hydraulic parameters of the overburden. This sensitivity analysis reduces the model ambiguity by elucidating the importance and dependence of such variables. The study outcomes aim to pave the way for the application of the approach in other postmining scenarios, whilst being an important step towards the simulation of more complex processes such as contaminant migration or groundwater rebound.

Geological setting
The Ibbenbüren coal-mining district is located in the northern part of North Rhine-Westphalia, north-west Germany (Fig. 2). Centuries-old mining industry was founded here on the exploitation of anthracitic coal seams encountered in a Carboniferous crustal block (Drozdzewski and Dölling 2018;Rinder et al. 2020). The coal-bearing sequence was brought to the surface by a compression event during the Upper Cretaceous, which produced a crustal uplifting of about 2 km compared to the Triassic and Jurassic foreland, currently under the Quaternary. (Drozdzewski 1985;Drozdzewski and Dölling 2018). This process created an island-like shape hill, with marginal faults limiting all its sides (Figs. 2 and 3). Compression forces also created antithetic faults that divided the block into two NNE-SSW striking horst structures separated by the Bockradenergraben. The latter represents a stratigraphic discontinuity of more than 250 m in relation to the rock sequence contained in the horst structures. (Lotze et al. 1962;Bässler 1970). For mining purposes, the two horst structures were named Eastfield and Westfield.
Lithology of the Westfield is well known based on geological documentation of mining activities. Rocks were deposited during the higher Westphalian C and lower Westphalian D stages, reaching a total thickness of more than 1,500 m (Blowes and Jambor 1990;Wüstefeld et al. 2017;Becker et al. 2019). The subhorizontal sequence follows an alternating pattern of sandstone and conglomerate layers, with sporadic presence of shales, wherein 81 coal seams are embedded (Lotze et al. 1962;Bässler 1970). This pattern is equally observed in the shallow overburden, where a fining-upwards array of sublithic sandstone layers (approx. 80%) appears interspersed with dark-gray shales (20%) and remnants of two coal seams (Bedoya-Gonzalez et al. 2021a). On the surface, thin Quaternary sediments and a few anthropogenic waste rock deposits occasionally cover the Carboniferous rocks.

Hydrogeological setting
Mining operations on the Westfield stopped in June 1979, with excavations as deep as 600 m below ground level (Klinger et al. 2019). After its closure, the area experienced controlled flooding up to 65 m above sea level (asl; Rudakov et al. 2014;Bedoya-Gonzalez et al. 2021b). At this elevation, groundwater reached an underground tunnel known as Dickenberg adit that runs along the southern edge of the field. Hydraulic connection from the Bockradener Graben to the adit is prevented with artificial underground dams. This generates the Dickenberg adit to drain the water from most of the Westfield, maintaining the phreatic level at 65 m asl. This includes the area enclosed by the Northern and Southern Carboniferous Marginal faults, Mieke Fault and Pommer-Esche Fault (Fig. 3), which represent effective hydrogeological boundaries (Rudakov et al. 2014;Coldewey et al. 2018). Furthermore, groundwater can be extracted only in areas higher than the Dickenberg adit, turning the terrain contour +65 m into an additional hydrogeological limit Since the phreatic level of the former coalfield is above the foreland surface (<55 m asl), percolated precipitation becomes the only groundwater input. The Geological Survey of North Rhine-Westphalia has calculated the groundwater recharge of the area as a percentage of the precipitation using a water budget model that considers seasonality, soil type, land use and vegetation (Herrmann et al. 2014). Outcomes show an intensive groundwater recharge in winter, with a maximum in December and January, where it can reach 40-80 mm/month. In November, February and March, the monthly sums of groundwater recharge are around 30-40% lower than the previous period. April and May represent transition months where significantly more water is evaporated due to the beginning of the vegetation period. During summer months, recharge drops to almost zero, reaching a maximum multiyear average of 10 mm/month. Finally significant recharge is observed again in October. In conclusion, the Westfield can be conceived as a well-defined closed system, where the discharge at the Dickenberg adit depends on the temporal distribution of the recharge and the percolation dynamics through the unsaturated fractured overburden.

The dual-continuum approach
This study implements TOUGH2 code (by Lawrence Berkeley National Laboratory) to formulate a dual-continuum model of the fractured and unsaturated overburden of the Ibbenbüren Westfield. The numerical solution of the software employs space discretization with integral finite differences and fully implicit first-order finite differences in time (Xu et al. 2000). Modeling scenarios include isothermal flow conditions, with air and water as the two phase components. The model requires the software to simultaneously solve two sets of flow equations at the same node: one for the porous medium and the other for the fracture network (Pruess and Narasimhan 1985). Here, the fractured system is represented as a porous medium by considering laminar flow through a network of small conduits. This assumption also accounts for the presence of undisturbed porous layers above the fractured medium as described for mining models.
If water inflow at the interface with the fractured continuum is greatly dependent on the hydraulic conductivity of the overlaying layers, the Reynolds number (Re) for such continuum (Eq. 1) would most likely be below 10, allowing turbulent flow to be neglected (Bear 1972): where ρ w is the water density (kg/m 3 ), K u-f is the hydraulic conductivity (m/s) at the interface between the undisturbed layer and fractured continuum, μ w the water viscosity (Pa s), and d is a representative length dimension for the porous medium (m), commonly taken as a mean grain diameter.
(1) Fig. 3 Geological map and cross-section showing the hydrogeological boundaries and Dickenberg adit. The A-A′ line in the map denotes the position of the cross-section, which is vertically exagger-ated 5× to better detail the shallow overburden structure (modified after Bedoya-Gonzalez et al. 2021a) Following this approach, the computational burden is reduced to solving a pair of Darcian flow equations rather than coupling, for example, a numerical Darcian solution to a Navier-Stokes equation (Chen et al. 2012;Aguilar-López et al. 2020). For an unsaturated system, flow can be described by two Richards equations (Richards 1931), coupled by a scaled exchange term (Γ ex ) that considers the interaction between both media (Gerke and van Genuchten 1993): where the subscripts m and f indicate matrix and fractured continua respectively (hereafter replaced by the subscript c for space saving reasons since the equations for each continuum look the same), θ is the specific volumetric water content (m 3 /m 3 ) defined as the product between the porosity (-) and water saturation (m 3 /m 3 ), h is the hydraulic head (m) and K is the hydraulic conductivity (m/s). K in turn can be obtained from Eq. (4) (Pruess et al. 2012): where g is the gravity acceleration (m/s 2 ), k c the absolute permeability (m 2 ) and k rw c a relative permeability coefficient (-) directly related with the water saturation of each medium. For saturated media k rw c is equal to 1, while under unsaturated conditions, the coefficient varies according to a number of factors described by the van Genuchten-Mualem parametric model (Mualem 1976;van Genuchten 1980): with S e c being the effective saturation (-), θ r c the residual volumetric water content (m 3 /m 3 ), α c a scaling fitting parameter (1/m), and n c and m c fitting parameters associated with the pore distribution. Finally, the exchange term (Γ ex ) in the Eqs.
(2) and (3) can be defined as (Gerke and van Genuchten 1993): and α * is the exchange coefficient (1/s), which depends on the size and shape of the mesh elements according to: In Eq. (8), β is a geometry grid factor (3 for rectangular blocks, 15 for spheres), a is the distance between the center of a matrix block and the adjacent fracture and γ w is a dimensionless scaling coefficient usually set to 0.4 (Pruess 1983;Gerke and van Genuchten 1993;Gerke et al. 2007;Kordilla et al. 2012)

Model calibration, efficiency and sensitivity analysis
The required parameters for fractured and matrix continua are manually calibrated using daily discharges recorded throughout the year 2008. The process is carried out through a trial-and-error scheme in which the simulated signal is visually adjusted with respect to the measured discharge. It starts with the adjustment of known parameters, including the permeabilities and porosities of the rocks, followed by the unsaturated properties of the porous units and hydraulic properties of the fractured continuum (usually reported in the literature and mining reports), and ending with the setting of the unsaturated parameters of the fractured continuum that have the lowest degree of knowledge. Once calibrated, the quality of the obtained parameters are validated using the measured discharges of the year 2017. Thereafter, the Nash-Sutcliffe efficiency (NSE) criterion is used to quantify the predictive power of the model for simulating the water discharge of the calibration and validation periods. This is a normalized coefficient that determines the relative magnitude of the residual variance compared to the measured data variance. NSE indicates how well the plot of observed versus simulated data fits the 1:1 line (Moriasi et al. 2015). While NSE equal to one represents a perfect match of the simulated to the observed data, values equal or lower than zero indicate that the model predictions are as accurate as the mean of the observed data. Thus, NSE is a convenient tool to evaluate the ability of a model to predict observations in long-term continuous simulations by comparing its effectiveness with respect to its mean value (Gupta et al. 2009).
A sensitivity analysis is additionally performed to determine those variables that can be set over a reasonable range, without affecting the efficiency of the model. The process is conducted by changing individual parameters of each continuum to at least four different points along the entire range of values. Although simple, this method appropriately quantifies the factors influencing the distributed regional model, maintaining narrow and reasonable control of the outcomes of interest with relatively few simulations (e.g., see Herman Rakovec et al. 2014;Devak and Dhanya 2017). Here, the lower and upper limits of each parameter are set on the maximum and minimum measured values present in the literature. The impact of the parameters is evaluated by calculating the root mean square error (RMSE) of the computed outflows compared to the calibrated discharge signal of the year 2008. The RMSE metric is preferred over the NSE in this step to overcome the known weaknesses of the latter in identifying differences in timing and magnitude of peak flows and shape of recession curves (e.g., see Krause et al. 2005;Pushpalatha et al. 2012;Moriasi et al. 2015). The calibrated signal of the year 2008 is also employed to facilitate interpretation regarding the best obtained adjustment, while preventing the computed RMSE values from being affected by outliers of measured daily discharges. Finally, interrelations between factors are evaluated with a codependency analysis, changing in pairs all the modeled parameters across the predetermined ranges.

Geometry and grid construction
The conceptualization of the modeled area is depicted in Fig. 4. Its design and dimensions are constructed according to the limits previously described in the hydrogeological section and illustrated in Fig. 3. The 3D system is discretized into vertical columns of 100 × 100 m 2 base and variable height to represent the concave topography of the Westfield. The columns' height decreases about 5 m per column according to the average thickness of the rock layers. This size considers actual measurements made on core samples of the shallow overburden in the study area (Bedoya-Gonzalez et al. 2021a). The total number of columns for each height interval is, then, normalized by intersecting the Fig. 4 Constructed grid for the DC Model. Recharge is applied directly onto the weathered blocks while an open boundary condition is set along the adit element. For each block colored in orange, the mesh is divided into two continua (fracture and porous media) with the same physical location digital elevation model (DEM) of the Westfield with the base boundary (i.e., 65 m asl). Finer horizontal discretizations of the columns are neglected since the model only computes gravity-driven vertical flow. Moreover, TOUGH2 requires fracturing information (i.e., spacing, number of sets, and shape of matrix blocks) and not block size to set the nodes location for calculating the driving pressure gradient at the matrix/fracture interface.
The average measured thickness of the layers is also used to subdivide the grid into blocks of 5-m vertical spacing (Fig. 4). Further vertical refinement is not implemented considering the rather homogeneous sequence, which is discretized into three lithologies, with predominance of sandstones (more than 80%). More layers would slightly modify the permeability of the columns due to minor saturation differences within the smaller blocks, but at higher computational cost. Like the thickness, horizontal layers are assigned within the columns maintaining the lithological proportions and stratigraphic position described in Bedoya-Gonzalez et al. (2021a). Exceptions are weathered intervals within the uppermost 5 m of each column and zones above 140 m, which correspond to anthropogenic rock deposits assigned to the weathered rock domain. Mining-induced fractures were estimated to extend up to 45 m above the uppermost mined seam, i.e. above the Dickenberg adit. This estimation is based on mathematical relationships between thickness of the excavated area and thickness of the overburden as previously established in Palchik (2003), Bai and Tu (2019), Guo et al. (2019) and Zha et al. (2020). However, additional grids with heights of 30 and 40 m are created to assess alternative values derived from these models. In any case, the base of the weathering zone turns into the uppermost limit of the fractured continuum. This causes the splitting not to extend to the first element of each column, even if the thickness of the overburden is less than the assumed height. Finally, the fracture density (i.e., number of vertical fractures per grid block) is assumed to be equal to 10, resulting in a total exchange area of 10,000 m 2 between both continua. As there is no study to the authors' knowledge that gives an approximation to this value, grid cells with 5 and 20 fractures per block are also constructed for the sensitivity analysis of the parameter.

Boundary and initial conditions
The lateral sides of each column are defined as no flow boundaries while the lower boundary is set to allow free drainage under gravity force. A specified flux boundary is set at the top of each column (i.e., in the weathered blocks) to account for diffuse recharge (Fig. 4). The amount of flux at the top is based on the recharge model developed by the Geological Survey of North Rhine-Westphalia (see section 'Hydrogeological setting'). To account for weather variability, daily precipitation measurements are normalized with respect to the precipitation values used in the recharge model. Precipitation above the coalfield is calculated using the arithmetic mean values of daily events measured at the two closer meteorological stations, Mettingen and Laggenbeck, for the years 2007 to 2017 (data provided by the LANUV NRW). As a result, recharge imposed on the model intrinsically considers the influence of the climatic and physical characteristics of the area, varying daily according to the precipitation events- Table S1 in the electronic supplementary material (ESM).
Under unsaturated conditions, TOUGH2 requires water saturation as the initial condition for each grid element. Here, saturation values are computed from a long-term simulation that starts with the residual water content of each lithology (see the following section). First, constant recharge of 240 mm/year, equivalent to the long-term average groundwater recharge in northwest Germany (BGR 2021), is assigned over 100 years (approximate time from mine opening to 2008). Daily recharge values for the years 2007 and 2008 are then applied to ensure the seasonal variability of the water storage during the model calibration. Finally, daily inputs are continuously applied for another 9 years to maintain the temporal variability for the model validation in 2017.

Parameterization
Parameters used in the model are listed in Table 1. These values were extensively searched in the literature and assigned according to the characteristics of the matrix and fractured continua (Freeze and Cherry 1979;Kordilla et al. 2012;Parajuli et al. 2017;Coldewey et al. 2018;Bedoya-Gonzalez et al. 2021a). For model calibration and sensitivity analysis, parameters are varied within reasonable data ranges consistent with known field conditions (values in parentheses Table 1). Since there are no documented values for the hydraulic properties at the interfaces, these are set equal to the harmonic mean between the two continua. This forces the local exchange to be dominated by the lower conductive domain (i.e., the matrix), which might be closer to reality (Ray et al. 2004;Song et al. 2018;Aguilar-López et al. 2020). Specific storage coefficients for both media are neglected since the model faces unsaturated conditions and water release by compression is irrelevant.
Given that parameters of the van Genuchten-Mualem function have been developed for porous media, values assigned to fractured continuum have questionable physical meaning (e.g., see Sauter 1992; Kordilla et al. 2012). Here, these parameters are allocated to recreate the expected hydrological response of a highly permeable medium. Total porosity is set to 0.99 for the software to assume open fractures without solid elements, as identified by Bedoya et al. 2021 in core samples of the study area. However, the influence of a solid phase within the fractured continuum (e.g., sediment mixing flow, clay smear or cement) is also considered within the sensitivity analysis by reducing the open porosity down to 50%. The residual water content (θ r ) is fixed at 0.05 to allow the continuum to transmit water during dry periods. If this value was zero, the relative permeability would nullify the exchange with the matrix, generating numerical insufficiencies. The small amount of water saturation acquires physical meaning, for example, if a process such as gravity-driven film flow on the fracture walls was to be considered (Tokunaga and Wan 1997;Kordilla et al. 2012). Finally, the parameters m and α also reproduce the behavior of a large-pored material (e.g., gravel) that will saturate and transmit water rapidly during rainfall events. However, α parameter needs to be adjusted during the calibration process.

Calibration and validation of the model
The parameterization within the model was calibrated using daily discharges recorded throughout 2008 and subsequently validated with data of 2017. Fluid flow was simulated in the vertical component of each column, allowing lateral exchanges between fractured and porous media. The total discharge of the study area was obtained by summing the water flux of the two continua at the lower boundary of every column (i.e., at 65 m asl). This is equivalent to a dual permeability model, in which the global flow occurs in both continua. In the case where columns are covered by waste rock deposits, discharge measurements were made around 140 m asl, considering that an auxiliary drainage directs the water from these deposits towards the Dickenberg adit (B. Nippert, RAG GmbH, personal communication, 2019). Figure 5a shows the measured and simulated discharge volumes of the Dickenberg adit for the year 2008. The graphed curve corresponds to the best match using the values listed outside the parentheses in Table 1. Model calibration was accomplished by fitting the observed and simulated discharge curves through a trial and error exercise that was performed in a series of more than 500 runs. Due to software limitations, it was not possible to perform, compare or supplement the calibration process with automatic inverse methods. Notwithstanding, manual calibration alone is expected to produce process-based (i.e., conceptually realistic) and reliable predictions from a clear understanding of the model structure (i.e., hydrological processes) and physical characteristics of the area (Ndiritu 2009;Arnold et al. 2012;Acero Triana et al. 2019). To prove it, the model reliability was validated for another 12-month period using the discharge of the year 2017 (Fig. 5b). Both, calibration and validation processes used variable time stepping. The initial time step was set at 100 s being systematically doubled if convergence occurred within the first four iterations up to a maximum of 1 day. Conversely, time step was reduced in multiples of 4 if the system exceeded 8 Newton-Raphson iterations to solve the connected equations. At the end the discharge was computed on a daily basis to compare the simulated signal with the measured data in the Dickenberg adit.
Modeling results also describe the temporal and spatial interaction between the fractured and matrix continua. For instance, the initial response of the system to a strong precipitation event in winter is a significant water transfer from the fractures to the matrix (Fig. 6a). However, much of this water returns to the fractures during spring and summer months (negative exchange values), as water pressure in the matrix increases, especially at interfaces with shale layers. On the other hand, Fig. 6b displays the exchange rate at the first and last fractured sandstone layer during a high rainfall period in March 2008. The delay of the inflow pulse to the matrix is associated to the pathway followed by the water within the fractured continuum and the storage capacity of the porous medium. Note that the exchanged water volume is much lower in the last layer because of its high saturation and spatial proximity to the water table. Besides, the exchange pulse shows a broader timeline in the last sandstone layer when compared to the first layer. This is considered a dispersive phenomenon which may be caused by the difference in hydraulic properties between the two continua leading to variations in percolation velocities across the overburden. Lastly, Fig. 7 shows saturation variations due to water exchange between fractured and matrix continua for one of the columns. Results reveal that the fractured continuum saturated and desaturated faster and in higher percentages than the matrix medium, although the amount of water either gained or lost in the latter is much higher than in the fractured one. Figure 7 also summarizes how the  . 6 Temporal and spatial distribution of the total fluid exchange and water saturation. a Total daily exchange between the two continua for the whole system (positive values represent exchange from the fractures to the matrix while negative values from the matrix to the fractures). b Spatial variation of the total water exchange between the first and last fractured sandstone layers matrix absorbs water from the fracture in the winter, when saturation is higher, returning a large amount in the summer season due to pressure increase, especially at the interfaces with the shale layers. Despite this, the closest matrix blocks to the surface continue to absorb water coming from the unfractured sequence above them, preventing the fractures from easily discharging percolating rainwater into the adit during the drier months.

Model reliability and sensitive parameters
The simulated signals for calibration and validation periods display good visual agreement with the measured data. Results show few discrepancies in the discharge peaks for the winter months and a small signal advance in time, likely derived from the nonestimation of the rapid flow from the adit to the measurement point. The reliability of the model in fitting the measured discharges is further proven by NSE values close to one (Fig. 8). The largest discrepancies from the 1:1 trend line correspond mostly to an overestimation of the simulated discharge (values above the trend line in Fig. 8) derived from the aforementioned signal advance. The steep rising discharge limbs in winter derive significant differences between the measured and simulated values. For example, there is a difference of about 6,000 m 3 between the two discharge signals at the end of January 2008, because the peak of the simulated signal is reached 1 day earlier compared to the measured discharge. The simulated signal also presents underestimations in the discharge calculation associated with recession limbs in winter. However, differences between the two models are much smaller than those observed with the rising limb and are associated to punctuated events. Table 2 provides an overview of all the field parameters required by the software as well as their sensitivity within the model. The range of values considered for each parameter matches those reported in Table 1, in addition to the variations of the fracture density and fractured zone height discussed in the section 'Geometry and grid construction'. Here, parameters have been categorized as sensitive if the maximum dispersion of the data points around the regression line of the calibrated mean annual discharge is greater than 10%, which is equivalent to an RMSE of 1,360 m 3 / day for the year 2008. As some parameters are likely to be insensitive when varied independently due to the complexity of the model, codependences are evaluated by varying two parameters at the same time. Thus, Table 2 also shows the highest RMSE discharge values for pairs of parameters with just one being sensitive (i.e., linear relationship) and pairs where both parameters are sensitive (i.e., nonlinear relationship). The analysis, however, neglects those recurrent insensitive and unrelated pairs of parameters associated Fig. 7 Change in water saturation and fluid exchange between the two continua at two temporal points. Water exchange is given in m 3 /day, while the variation in water saturation is measured with respect to the initial water content at the time the daily variable recharge is applied (i.e., the water saturation at the end of the 100-year long-term simulation) with the model setup. For example, the influence of shale permeability when varied with respect to the porosity of the fractured continuum is considered redundant for two reasons: (1) fracture porosity minimally influences the already affected discharge signal by the fracture permeability (i.e., there is a linear relationship between the two parameters where only the fracture permeability influences the discharge) and (2) fluid exchange at interfaces was set up to depend largely on matrix permeability instead of fractured continuum permeability (see section 'Parameterization'). As a result, the behavior and RMSE values are almost identical to those already reported for scenarios in which the shale permeability is varied along with the permeability of the fractured continuum.
Complementarily to Table 2, Fig. 9 displays the discharge curves of individual parameters with the highest sensitivity, while Fig. 10 shows the types of codependences that exist between some couples (all the codependences between parameters can be found in Fig. S1 of the ESM). From Fig. 9, it is observed that shorter water-conducting fractured zones (Fig. 9a) and higher sandstone permeabilities (Fig. 9b) broaden and enlarge the recession limbs, with up to 25% more water discharge. On the one hand, shorter fractures generate longer water residence times in the uppermost pristine layers, while higher sandstone permeability enhances the exchange from the fractured to the porous medium. Note that the exchange process between the two media is dominated by the permeability of the matrix (see 'Parameterization'). This configuration causes strong discharge peaks in winter months when the sandstone permeability is reduced (Fig. 9b), as well as minimal signal differences after increasing the fracture permeability (Fig. 9c). The insensitivity of the model to higher fracture permeability is also related to the placement of a weathered layer on top of the columns, which acts as a buffer for the amount of precipitation that infiltrates into the continuum. On the other hand, lower fracture permeability results in higher base flow discharges at the expense of decreased winter peaks (Fig. 9c). Observe that for the lowest value (i.e., k f = 1 × 10 -12 m 2 ), the fractured medium reaches a permeability similar to that of the porous medium, generating the closest scenario to an EPM model. Finally, it was observed that the alpha (α) parameter of the van Genuchten-Mualem function is directly related to the magnitude of the discharge (Fig. 9d). Higher values cause the fractured continuum to retain less water during the winter months, transmitting it more quickly and generating greater discharge.
The analysis of codependences revealed 14 combinations in which both parameters were shown to be sensitive. For example, the simultaneous variation of sandstone permeability (k s ) and porosity (Φ s ) shows pronounced fluctuations, especially for low values (Fig. 10a). Whereas for the calibrated permeability (i.e., k s = 5 × 10 -14 m 2 ) Φ s is insensitive; a k s of one order of magnitude lower turns the porosity sensitive with a RMSE of around 2,500 m 3 /day for 8% and 1,400 m 3 /day for 15%. A more abrupt change is observed when varying the fracture alpha (α f ) parameter over the lower limit of the k s (Fig. 10b). The lowest values of both parameters produce a minor sensitive combination, with a RMSE of about 1,500 m 3 /day. In contrast, the combination of the highest α f with the lowest k s increases the RMSE to almost 5,000 m 3 /day. It was also observed that the α f exhibits mostly a linear relationship with almost all parameters evaluated. For example, Fig. 10c shows how the RMSE of the discharge varies widely as α f increases without being affected by the residual water saturation of the fracture (θ r f ). Similar linear relationships are observed for couples that include k s and k f . Lastly, Fig. 10d illustrates how Φ s and θ r s are insensitive over the entire evaluated spectrum. Other recurrent insensitive variables are the van Genuchten m parameter of the porous layers, the permeability of the shales, the fracture continuum porosity (Φ f ), and the fracture density. The behavior of the latter is surprising as higher densities generate larger exchange area between the two continua. The low permeability of the matrix seems to be a highly restrictive parameter in this model, maintaining the global discharge of the study area. This same trend would be partially responsible for the low sensitivity of the fracture continuum porosity. The reduction of Φ f leads to higher water saturation in the continuum, increasing its potential to the matrix. However, matrix tightness causes minimal variations in the amount of exchangeable water for any permeability assigned

(d)
to the sandstone layers (see Table 2). On the other hand, the increase in water saturation is not sufficient to increase the water flow and discharge within the fractured continuum. The small saturation changes are a consequence of the storage capacity and low transmissivity of the pristine rock layers above the fractured zone. As a result, the van Genuchten-Mualem parameters employed for the fractured continuum derives minimal variations of the relative permeability values. Fracture porosity would, however, play a more significant role in a reactive model in which the solid fraction interacts chemically with the percolating water.

Discussion
The good agreement between the measured and simulated discharge of the Ibbenbüren Westfield suggests that the dual-continuum approach is appropriate to characterize water movement through the heterogeneous and disrupted overburden of underground mines. The constructed model reproduces the bimodal behavior of the adit discharge taking into account the rapid flow of a fractured continuum, the storage capacity of a matrix continuum and the water transfer between both media. It was not possible to model a comparable signal trend or magnitude by assigning similar hydraulic properties to both continua under the same water balance or overburden arrangement. The dashed lines in Fig. 9b,c show how the model averages the discharge of the adit when the permeabilities of both media are close to each other (i.e., the closest scenario to an EPM model). Neither the winter discharge peaks (discharge is underestimated) nor the recession during the summer months (discharge is overestimated) are correctly simulated, resulting in scenarios with the largest discrepancies compared to the measured data. This behavior could be one of the reasons why the average discharge of the adit was used to calibrate the most recent EPM model of the Westfield (Coldewey et al. 2018). In contrast, the simulated signal where the permeability of the fractured medium is significantly higher than the porous one (i.e., solid lines of the Fig. 9b,c) would support the convenience of the DC approach over an EPM. The importance of the fractured medium for the hydrology of the Westfield had, in fact, been previously identified in mining reports. For instance, Lotze et al. (1962) and Bässler (1970) observed a direct correlation between seasonal precipitation and the amount of discharged water for works deeper than 200 m. Both authors also described a sharp drop in the water table for the whole field when mining was active. The development of a continuous cone of depression suggests a direct connection between the surface and underground works considering the low permeability of the sandstones, usually below 1 × 10 -12 m 2 (Bässler 1970;Wüstefeld et al. 2017;Becker et al. 2017). More recently, Bedoya-Gonzalez et al. (2021a) studied several drill core intervals segmented by subvertical fractures just above the Dickenberg coal seam. The presence of weathering fronts at both sides of the fracture planes together with the absence of healing elements suggest a recent origin probably linked with mining activities.
Simulations showed that the permeable fractured continuum dominates the flow system in the winter months, facilitating rapid responses after precipitation events. This behavior generates for example the two strong discharge peaks in February and March 2008 (Fig. 5a). For both events, the volume and subsequent recession limb were effectively simulated by the model, producing a good match with the measured data. The greatest model mismatch occurred in the length of the first discharge peak, which could not be fitted with the calibrated parameters. Initially, this inconsistency may be related to the type of precipitation. A similar factor could cause the slight discrepancies observed in the small peaks of January 2008 and 2017. It is worth noting that only liquid precipitation was considered when applying the daily recharge to the model. In the case of strong snowfall, recharge will be conditioned to the snowmelt over the following days and weeks, generating a delay in the discharge as well as long-lasting peaks during the melting period (Buttle 1989;Flerchinger et al. 1992;Eckhardt and Ulbrich 2003). If so, the recharge adopted for February 2008 was likely underestimated. A more detailed local water balance would be necessary to quantify this effect, as it is clear that the system reacts easily to recharge events due to the fractured continuum. Other reasons for the aforementioned discrepancies may include variations in the layer geometry, layer thickness or a more heterogeneous overburden, all of which were excluded during the model calibration.
Diversely, the regional recession limb of the discharge signal is influenced by the matrix flow as shown in Fig. 6a. During summer months the saturation of the fractured medium drops, increasing the relative potential of the matrix. The difference causes the water contained within the porous material to be slowly released back to the fractures and towards the adit. The continued desaturation of the system during these months increases the capillary pressure of the matrix sucking in fresh precipitation water. As a result, the system does not readily respond to precipitation events in summer. This phenomenon may be also affected by the presence of large pores. Be reminded that this model approximates the fractured continuum as a porous medium with enhanced hydraulic properties, similar to a gravel bed. Under low saturation conditions, large pores behave more like a barrier where air cannot be easily displaced by new percolating water (Villagra-Mendoza and Horn 2018; Aguilar-López et al. 2020). Thus, the employment of the Richards equation for both continua is a common simplification that fits the study quite well when compared to other DC models. Simulations of karst systems, for example, often adjust boundary conditions to include the nonlaminar fast component of conduits (e.g., Kordilla et al. 2012;Dal Soglio et al. 2020). This issue seems to be related to the rapid percolation process of meteoric water through wide and continuous karst elements reaching the surface. Within coal mining areas, the occurrence of an undisturbed overburden above the fracture zone reduces the flow velocity of water entering into the fractures, favoring the application of a Darcian porous-media numerical solution for the entire system. Values from 0.0028 to 0.28 corroborate the laminar flow regime when calculating the Reynolds numbers (Eq. 1) at the upper interface of the fractured continuum with the sandstone layers and using their highest hydraulic conductivities (i.e., fully saturated conditions). Complementarily, the development of shorter, narrower and rougher walled fractures allows the fractured continuum to more closely resemble a porous medium, as it is affected by matrix potential and capillary forces (Liu and Bodvarsson 2001;Finsterle et al. 2003). This is added to the fact that change in fracture density along the vertical direction would lead to only part of the elements becoming active to water flow, increasing tortuosity and offsetting the potential fast component (Liu et al. 1998(Liu et al. , 2003. The sensitivity analysis proved major variation of the discharge curve when varying the permeability and van Genuchten-Mualem parameters of both continua; however, the model ambiguity was remarkably reduced by having a strong a priori knowledge of the parameters for the matrix continuum. The hydraulic properties of the sandstones were obtained from local investigations conducted at the Ibbenbüren Westfield. The mining framework was an advantage for the model setup since mining companies typically hold considerable amounts of geological information collected during the exploration and exploitation stages. There are other sensitive parameters though, with scarce reported values or minimal information about their measurement. One of them is the α parameter of the van Genuchten-Mualem function for the fractured continuum. Although many studies did use it to simulate unsaturated flow in DC and discrete models, only a few have explained its meaning for fracture networks and respected uncertainty for nonporous materials (e.g., Vogel et al. 2000;Liu et al. 2003;Colombo et al. 2013;Koohbor et al. 2020). In this study, all the parameters of the van Genuchten-Mualem function were initially assigned to correspond to a permeable medium with large pore diameters; however, the α value had to be adjusted to a value similar of the sandstone layers to bring the system and adit discharge to a reasonable state. The adjustment was done after constraining many of the matrix parameters with real data. It would be useful to determine whether this calibrated α value has a physical meaning, i.e. whether or not it could be adapted for other model applications, or it is just a fitting factor. Fractured zone heights represent another sensitive feature that must be carefully set. The parameter value varies from site to site, depending on technical and geological factors such as mining depth, lithology, layer thickness and rock deformation (Meng et al. 2016). A more accurate characterization of each postmining overburden would help to improve the model outcomes. Lastly, even though fracture density was highlighted as an insensitive variable, its influence should not be neglected, especially when working with more permeable rock sequences or contaminant transport. Here, the number of fractures and their aperture were calculated using the cubic law and calibrated permeability. For example, if a grid cell of 1 m 3 contains a single set of fractures with permeability of 1 × 10 -10 m 2 , the fractured continuum can be either depicted with 10 elements of 1.1 × 10 -5 m aperture or 20 elements of 7.7 × 10 -6 m aperture. The difference, then, lies in the contact surface between the two continua. For the second case the surface area is doubled compared to the first. This generates a significant increase in the volume of rock exposed to water and air due to fractures, which would greatly influence the generation and transport of contaminants. The inclusion of the latter would also help to evaluate the dispersive capability of the model, complementing what was observed in the infiltration pulses in Fig. 6b. The concept of hydrodynamic dispersion treated as a Fickian diffusion process is diminished in TOUGH2. However, a similar and more realistic process could be obtained by considering molecular diffusion together with flux dispersion due to mesh heterogeneities. The latter was apparently achieved in this model after considering various hydraulic and physical properties for each continuum.

Conclusions
The application of a dual-continuum model successfully reproduces the transient discharge of the Dickenberg adit in the Ibbenbüren Westfield. The good agreement between the calculated and measured signals for the years 2008 and 2017 was possible by coupling a high permeability fractured continuum with a low-conductivity matrix continuum. The unified water flow through the unsaturated overburden was simulated using two Richards equations, linked together with a scaled exchange term. Results show the fractured continuum as the main factor responsible for strong and short discharge events during the months of heaviest precipitation, while flow within the matrix continuum greatly influences the smooth recession limb of the dry period. The direction and quantity of exchanged water between them follow pressure differences according to the dominant flow element in each season (i.e., mainly from the fracture to the matrix during winter months). The largest discrepancies in the results were associated with the regional recharge model used as input to the simulations. The lack of a more precise water balance that additionally includes, for example, the type of precipitation, altered the shape of some discharge peaks in winter due to the rapid reaction of the fractured medium.
Results of the sensitivity analysis for the Ibbenbüren Westfield showed the permeability and van Genuchten-Mualem α parameter as two of the most sensitive variables. However, a solid a priori knowledge of the study zone helped to restrict the ranges of porous media values, reducing the ambiguity of the model. This allowed fitting the fracture continuum α parameter during the calibration process, for which sparse data were available. The same served to fix the fractured zone height and fracture density variables. The construction and application of a similar adjustment process for others DC models (i.e., different coal mines) would be required to know whether the sensitive variables fitted here remain valid or change following variations in hydraulic parameters, layer heterogeneity or external recharge factors. In this sense, the existence of extensive databases in many coal mining areas would allow for relatively easy construction, calibration and confirmation processes. However, the acquisition of more site-specific data and the inclusion of the transport component would still be essential to reduce the uncertainty of the fractured continuum parameters. The transport component would in fact help to extend the analysis of the dispersion phenomenon obtained from the heterogeneity of the two continua when assessing the production and migration of contaminants. In contrast, the simulations of the study area showed little sensitivity to changing the porosity of the matrix continuum, as well as the m and n van Genuchten-Mualem parameters for both continua. Therefore, these variables can be fixed by reasonable estimates.
Funding Open access funding provided by Paris Lodron University of Salzburg. This work was financially supported by the Forum Bergbau und Wasser-FBW.

Conflict of interest The authors declare no conflict of interest
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/.