A Multiple Interactive Continua Model (MINC) to Simulate Reactive Mass Transport in a Post-Mining Coal Zone: A Case Study of the Ibbenbüren Westfield

We tested the suitability of the multiple interactive continua approach (MINC) to simulate reactive mass transport in a disturbed post-mining coal zone. To the authors’ knowledge, this approach has not been employed in such mining settings despite its relative success in other environmental fields. To this end, TOUGHREACT software was used to set up a MINC model of the unsaturated overburden of the Ibbenbüren Westfield. With it, we examined and evaluated water–rock interactions in both the fractured and porous continua as the main driver of elevated hydrogen, iron, sulfate, and chloride concentrations in the coal mine groundwater. Long and seasonal geochemical signatures were obtained by formulating and applying a five-stage modelling process that depicts the mining history of the area. The simulation results agree well with the concentrations and discharge trends measured in the mine drainage. Oxygen and meteoric water flow through the fractured continuum, leading to a high and steady release of hydrogen, iron, and sulfate ions derived from pyrite oxidation in the matrix continua closest to the fractures. Likewise, high chloride concentrations resulted from the mixing and gradual release of relatively immobile solutes in the matrix as they interacted with percolating water in the fracture. In both cases, the use of a multiple continua approach was essential to resolve sharp gradients for advection and faster kinetic reactions, while reducing the model’s dependence on block size for diffusive transport at the fracture–matrix interface. The model further allows for the calculation and analysis of solute exchange and transport in the unsaturated overburden resulting from rebound and imbibition processes, something pioneering when compared to other models in the field.


Introduction
Environmental impacts generated by poor water quality from coal mining areas are well documented around the world (Banks et al. 1997;Doumas et al. 2018;Nieto et al. 2007;Wolkersdorfer et al. 2022a). Drainage of saline water contained in rock formations increases the concentration of chloride in freshwater bodies, causing quality detriment and toxicity to their communities (Elphick et al. 2011;Timpano et al. 2015;Turek 2004;Vöröš et al. 2021). Weathering of naturally occurring sulfide ores leads to the formation of acid mine drainage (AMD), which results in an abundance of iron and sulfate ions and the mobilization of potentially toxic metals (Acharya and Kharel 2020;Cravotta and Brady 2015;Gombert et al. 2018;Zhang et al. 2021). The eventual precipitation of hydrous ferric oxides tend to clog riverbeds, generating negative effects on the flora and fauna (MacCausland and McTammany 2007;McKnight and Feder 1984).
The drainage of the former anthracite coal mine in the Ibbenbüren Westfield is no exception to the aforementioned contaminants. Pumped mine waters from the field are acidic and characterized by high iron, sulfate, and chloride content, which often surpasses the harmful levels for freshwater ecosystems (Bässler 1970;Rinder et al. 2020). This situation has generated public interest in the long-term liabilities related to the treatment of the mineralized water since closure of the mine in 1979. There is a growing emphasis on both mitigating negative environmental impacts and recovering/reusing resources that result from mine drainage (Moodley et al. 2018;Simate and Ndlovu 2014;Wolkersdorfer et al. 2022b). The latter is becoming increasingly important due to the European Green Deal's focus on recycling as an important factor of the future economy. Either way, choosing the best treatment option or deciding on the possibilities of reusing the solid or aquatic resources requires a long-term prognosis based on the correct characterization of the elements of the hydrogeological system. However, the source and evolution of the contamination at Westfield has been widely debated but only partially characterized. Much of the attention has focused on identifying and dissolving efflorescent salts (i.e. solid by-products of pyrite oxidation) that may remain in empty spaces such as tunnels and mine faces. The work is done through box models, with little concern for the water-rock interactions that lead to their formation (e.g. see Klinger et al. 2019). Moreover, the approach is not feasible for longstanding generation and elevated concentrations of dissolved iron and sulfate after mine flooding, as occurs in the study area. Alternatively, Bässler (1970) attributed the high ion levels to the continuous and homogeneous interaction of meteoric water with the solutes and minerals contained in the bedrocks. In such a case, the characterization of the geochemical interactions could resemble an equivalent porous medium model, commonly applied in coal mining areas. The latter represents an oversimplification of the medium since the moderate permeability of the coal overburden limits the flow and interaction of groundwater with the matrix components. In reality, the water-rock interaction would be more complex, being closely linked to the creation of fractured zones derived from mining.
The importance of considering both fractures and rock matrix has been corroborated by Bedoya-Gonzalez et al. (2022) who used a dual continuum model (DCM) to fit the transient mine adit discharge of the Ibbenbüren Westfield. Their results indicate that the fractured continuum is primarily responsible for the intense and short discharge events throughout the months of heavy rainfall. Meanwhile, the flow in the matrix continuum significantly affects the gradual recession limb of the discharge during dry periods. Regarding the geochemical component, development of conductive fracture networks results in new and extensive exposure of rock components to atmospheric conditions for weathering. This process has been reported to be much more pronounced in shallow aquifers above mined panels, where groundwater drawdown enhances chemical reactions with atmospheric oxygen and percolating meteoric water (Booth 2002;ERMITE-Consortium et al. 2004). Chemical reactions typically include dissolution of carbonate minerals, reductive dissolution of metal hydroxides, and oxidation of sulfide minerals, which mobilize cations, anions, and metals from rock strata (Jankowski and Spies 2007). Moreover, fractures serve as the primary flow pathways, with advection being the main transport mechanism, whereas transport within the adjacent rock matrix is notably slower, mainly driven by diffusion. The latter may lead to movement and mixing of preexisting waters along the rock sequence. The effect of this mixing, either a decline or an enhancement in water quality, hinges on the chemistry of the recharged source (Booth and Bertsch 1999;ERMITE-Consortium et al. 2004).
The objective of this study was to evaluate the suitability of the multiple interactive continua (MINC) approach for the characterization of multispecies mass transport and a range of geochemical reactions in geomechanically disrupted coal mining areas. Typically, intensive coal mining generates a high density of fractures in the overburden that would allow them to be grouped within a continuum, assigning it common properties distinct from the surrounding porous units. Grouping all fractures as a separate continuum represents a significant advantage of the MINC approach over equivalent porous media models, which usually fail to represent the dual characteristic behaviour of fractured systems. Moreover, the MINC approach requires less data than a discrete fracture model, making it more practical. To test the approach, this study employed TOUGHREACT software to set up a 2-D cross-sectional MINC model of the unsaturated and partially fractured Ibbenbüren Westfield overburden. The construction process extended the already calibrated DC model of the flow component to a model with five strings of nested continua: one for the fractures and four for the rock matrix. The further subdivision of the matrix enabled us to describe chemical and multispecies mass transport gradients within the matrix blocks during their long transient periods of interporosity flow (weeks and months). For the model assessment, this study simulated and compared the longterm and short-term contaminant formation, deposition, and transport in the study area since the start of mining operations. For this purpose, dissolution, precipitation, and dilution processes related to iron, sulfate, and chloride-bearing minerals and solutions were modelled using a five-step process that emulates the mining history of the area.

Geological and Hydrogeological Context
The study area is part of the Ibbenbüren mining district, located in northwestern Germany (Fig. 1). It is one of the two collieries that mined anthracitic coal from a Carboniferous crustal block. The latter was uplifted to the surface by the inversion of the Osning fault system during the Late Cretaceous (Drozdzewski and Dölling 2018;Rinder et al. 2020). As a result, the Carboniferous sequence became separated from the surrounding Triassic and Jurassic units (now beneath Quaternary deposits) by offsets of up to 2000 m (Fig. 2). Compressional forces further produced antithetic faults that broke the Carboniferous block in two horst structures (Dickenberg and Schafberg hills) and one SSW-NNE graben (Fig. 2). The Schafberg and Dickenberg hills were respectively labeled as East-and Westfield for mining purposes.
The lithology of the Westfield is well-documented through geological records of mining activities. The Carboniferous rock sequence reaches thicknesses of more than 1500 m, but mining only occurred at depths down to −170 m above sea level (a.s.l.) (Bässler 1970;Becker et al. 2019;Rinder et al. 2020;Wüstefeld et al. 2017). Within this interval, sublithic sandstone and conglomerate layers make up approximately 80% of the strata, while the remaining 20% is composed of shales and sporadic coal seams (Bedoya-Gonzalez et al. 2021a;Lotze et al. 1962). Occasionally, the rocks are covered by narrow beds of Quaternary sediments and anthropogenic waste rock deposits of the mined sandstones and shales.
During the Westfield operation, mine dewatering was largely influenced by the seasonal rainfall pattern. About 90% of the extracted water (10 m 3 /min) agreed with the volume of infiltrated precipitation. Such a great amount was possible by combining the absence of a Quaternary cover with the occurrence of anthropogenic water-conducting fracture zones. The remaining water volume was hypothesized as inputs from the surrounding Mesozoic aquifers and formation water (Klinger et al. 2019). The latter would represent up to 2/3 of the sandstone porosity volume, which ranges from 10 to 15% (Becker et al. 2017). Deep percolation of rainfall was also confirmed by the chemistry of the pumped water. Samples extracted along the 200 m working depths in the Westfield were classified as Ca-Mg-Na-SO 4 type (Bässler 1970). This composition is markedly different from the Na-Cl water type identified in seams at −400 m.a.s.l. in the Eastfield . However, pumped waters at the Westfield had pH values around 3 and unusually high concentrations of iron (100-1000 mg/L), sulfate (1000-3000 mg/L), and chloride (up to 200 mg/L) whose origin is not clear yet. In 1979, the Westfield colliery ceased production, leading to a gradual shutdown of the dewatering system. This caused groundwater to rebound over the next four years, reaching a level of 65 m.a.s.l (i.e. an average depth of 60 m below the ground surface). At this depth, the Dickenberg adit passively collects all of the mine water, directing it to a sequence of settling ponds for treatment (Bedoya-Gonzalez et al. 2021b;Rudakov et al. 2014). Initial geochemical signatures of the rebounded water in 1982 yielded discharges with iron concentrations up to 2000 mg/L, sulfate up to 4500 mg/L, and chloride up to 180 mg/L (Klinger et al. 2019). Six years later, iron and sulfate concentrations were significantly reduced to 500 and 2800 mg/L respectively, while chloride remained constant. From 1988 to 2018, the concentration of the three ions slowly continued to decrease, reaching annual averages of 1800 mg/L for sulfate and 180 mg/L for iron and chloride. Since the Dickenberg adit has imposed a phreatic level that is above the surrounding ground surface (i.e. < 55 m.a.s.l.), the quality and quantity of discharged water after the rebound depend on the rainfall temporal distribution and its seepage paths within the unsaturated fractured overburden. The above applies to areas inside the Southern and Northern Carboniferous marginal faults, Mieke Fault, and Pommer-Esche Fault (see Fig. 2), all of which are considered to be effective hydrogeological boundaries (Coldewey et al. 2018;Rudakov et al. 2014).

Fig. 2
Map and cross-section of the Westfield displaying the geology, hydrogeological boundaries, and position of the Dickenberg adit. The cross-section has been exaggerated vertically five times to better illus-trate the overburden, and its location is marked by the A-A′ line on the map (from Bedoya-Gonzalez et al. 2022)

Geochemical Processes
Of particular concern at the Westfield is the production of AMD, with a high load of iron and sulfate ions. Consistent with the mineralogical composition of the respective rocks (see Bedoya-Gonzalez et al. 2021a), groundwater composition would be dominated by the oxidation of small amounts of pyrite embedded within the rock matrix. The reaction progresses into the sandstone layers through the development of permeable sub-vertical fracture networks. The horizontal exchange of fluids from the fracture into the matrix allows meteoric (oxygenated) water to interact with pyrite, oxidizing it: For every mole of pyrite, reaction 1 yields one mole of Fe 2+ (aq) , two moles of SO 2− 4(aq) and two moles of H + (aq) (Nordstrom 1982). The oxidation process can alternatively occur through the following reaction: The control of the respective reaction pathway on the pyrite dissolution kinetics depends amongst others on pH, with Eq. (2) dominating under acidic conditions (Singer and Stumm 1970). Furthermore, reaction 2 is faster than reaction 1. Rate constants of pyrite oxidation by O 2 range from 1 × 10 −9 to 5 × 10 -10 mol m −2 s −1 , while those by Fe 3+ (aq) range from 1 × 10 −8 to 9.6 × 10 -9 mol m −2 s −1 (Nordstrom and Alpers 1997;Williamson and Rimstidt 1994). Nevertheless, these faster values result from the catalytic effect of ironoxidizing bacteria such as Acidithiobacillus ferrooxidans, a process not addressed in the simulations of this study. Under abiotic conditions, the oxidation rate of Fe 2+ (aq) to Fe 3+ (aq) happens on the order of 1.5 × 10 -10 mol m −2 s −1 , which makes it the limiting step of the reaction and thus yields similar pyrite oxidation rates to those reported with oxygen (Singer and Stumm 1970;Xu et al. 2000). This constraint is added to the fact that meteoric water could not carry high loads of Fe 3+ (aq) because of its instability at circumneutral pH, as is the case for the rainwater in the Westfield (Pierini et al. 2002). Once pyrite is oxidized, Fe 2+ (aq) can be transformed to Fe 3+ (aq) (reaction 3) to subsequently precipitate as amorphous iron hydroxide (reaction 4) within minutes (Hem and Cropper 1962). (1) The progress of reactions 3 and 4 in deep unsaturated zones would be constrained by the available oxygen remaining after the oxidation of pyrite. This can produce high amounts of Fe 2+ (aq) , SO 2− 4(aq) , and H + (aq) ions in solution within the porous medium. Alternatively, iron hydroxides can be converted to either goethite, which is thermodynamically stable at pH 2.5-5.5, or ferrihydrite/schwertmannite, stable at pH > 5.5, over aging periods involving months (Knorr and Blodau 2007;Van Der Woud and De Bruyn 1983). In the absence of carbonates in the area, the only probable mineral that can serve as a pH buffer is the kaolinite present within the rock matrix (Bain et al. 2001;Spiessl et al. 2007). Under acidic conditions, kaolinite dissolves according to the following reaction (Wieland and Stumm 1992): There are some documented cases of mine drainage waters exhibiting high levels of dissolved silica and aluminium (e.g. Blowes et al. 1992;Glynn and Brown 1996;Morin et al. 1988). Accordingly, this study also evaluated the pH buffering capacity of the kaolinite present within the rock matrix.

Model Setup
We used the multiple interacting continua approach (MINC Pruess 1983;Pruess and Narasimhan 1985) integrated into the TOUGHREACT software to set up a reactive transport model for the unsaturated overburden of the Ibbenbüren Westfield (For a detailed explanation of the mathematical principles underlaying the MINC modelling employed in this paper, please refer to the Supplementary Information). The model was constructed using the dual continuum fluid flow model developed by Bedoya-Gonzalez et al. (2022) as a basis. Specifically, all the calibrated hydraulic parameters defined in the original dual model were maintained, as well as the model geometry. The MINC model, however, divides the matrix blocks of which the mesh geometry is composed into four nodes, assigning to each of them different pressure, temperature, saturation, and concentration values. While all fractures were grouped into continuum #1 based on well-defined physical and chemical properties, continuum #2 includes the matrix material within a specific proximity to the fractures, and continua #3 and 4 comprise the matrix situated at successively further distances (Fig. 3). This subdivision allows the model to account for kinetic and diffusion reaction fronts among fractures and matrix. For the model's geochemical component, the study relies on the detailed mineralogical information reported in Bedoya-Gonzalez et al. (2021a). This mineral assemblage was allowed to interact either kinetically or at equilibrium with the irreducible formation water, percolating meteoric water, and rebound water over variable timespans. Furthermore, dilution and diffusion processes resulting from the interaction between waters and atmospheric oxygen were considered. The following subsections specify the features considered in the model construction.

Geometry and Grid Construction
The grid design and dimensions follow the marginal faults and well-defined boundary conditions previously explained in the hydrogeological section (Fig. 2). To account for the concave topography of the Westfield, the system is partitioned into vertical columns with a base area of 10 × 10 m 2 and variable height. The height of each column was reduced by 5 m based on the average thickness of rock layers measured on core samples of the shallow overburden (Bedoya-Gonzalez et al. 2021a). The same thickness was applied to divide the mesh into slabs with 5 m vertical spacing (Fig. 3). Only one column was constructed per terrain height due to the symmetry of the area and the fact that vertical infiltration was the sole component of the global flow towards the adit. Lithological and stratigraphic arrangement of the horizontal layers inside the columns was done following the position and proportions of the bedrock strata outlined in Bedoya-Gonzalez et al. (2021a). Anthropogenic deposits located above 140 m a.s.l. were considered to be weathered layers (Fig. 3). We assumed that mining-induced fractures within the overburden can reach a height of 45 m over the upper mined coal seam (i.e. above the Dickenberg adit). Fracture vertical growth follows mathematical relationships based on the composition and depth of the overburden, together with the thickness of the mined coal seam (Bedoya-Gonzalez et al. 2021a). However, fractures do not extend into the first element of each column (i.e. upper 5 m), which was assigned to a weathered layer. The resulting fractured overburden was discretized into one fracture continuum and four nested matrix blocks: the first three (closest to the fracture) are fixed slabs with thicknesses of 0.1, 0.4, and 1 m, while the fourth accommodates the remaining length to the next fracture. Vertical fractures were placed equidistant, with a linear density of 1 fracture per 10 m-long block, and a total interface area of 100 m 2 between continua. This horizontal 1-D array of nested grid blocks enables gradients for solute flow and transport within the matrix, perpendicular to the fracture.

Hydrodynamic Parameters
Different sets of permeability, porosity, and van Genuchten parameters were assigned to the sandstone, shale, and weathered lithologies according to their actual characteristics on the field (Table 1). Regarding the fractured continuum, values were assigned to reproduce the expected hydraulic behaviour of a high-permeability medium (e.g. gravel). These values emulate the real features of post-mining fracture zones, while meeting the software requirements for Constructed grid for the shallow overburden of the Ibbenbüren Westfield. The weathered blocks receive direct recharge, while each orangecolored block is partitioned into five continua that share the same physical location: one for fracture (Fr) and four for the matrix media (M n ) its modelling. For instance, the fractured continuum was given a porosity of 0.90 to enable the software simulating open fractures, such those discovered in the drill core of the study area by Bedoya-Gonzalez et al. (2021a). Moreover, to ensure water transmission during dry periods, the residual water content (θ r ) was maintained at a fixed value of 0.05. Setting θ r to zero would result in the relative permeability nullifying the exchange with the matrix, leading to numerical deficiencies. Note that the parameters used in both the matrix and fracture continua correspond to calibrated values for the best volumetric discharge agreement with respect to the measured data for the Dickenberg adit (Bedoya-Gonzalez et al. 2022). Therefore, disturbances of the flow component derived from porosity and permeability changes by dissolution and precipitation of minerals were not considered in this reactive model.

Boundary Conditions
Each column was assigned a no-flow boundary condition to their left, right, front, and back sides, while open boundaries were set on the interfaces between the fractured and four matrix continua. To account for diffusive recharge, a flux boundary was established in the weathered blocks located at the top of each column. The amount of water imposed was assumed to be either constant or daily variable depending on the modelling phase. A constant recharge of 240 mm/ yr was assumed to calculate a long-term water equilibrium. This value corresponds to the annual groundwater recharge average in NW Germany (BGR 2021). Alternatively, daily precipitation values were used to evaluate weather variability using the recharge model of the North Rhine-Westphalia Geological Survey (Herrmann et al. 2014). Above the weathered blocks, an atmospheric boundary with a constant oxygen pressure of 0.2 bar was imposed, enabling its diffusion and advection in the gas phase.
In most of the simulations, the lower boundary (i.e. Dickenberg adit) was set to permit free gravity drainage. However, the boundary was changed for a short period to a fully saturated Dirichlet-type to allow capillary equilibrium with the unsaturated sequence during groundwater rebound. It is important to point out that the simulated discharge and solute concentration were calculated at this boundary by computing the fluid and solute flux of the five continua with the adit. These values were normalized for the whole area by multiplying the flux of each column by its areal representation according to the digital elevation model of the Westfield. As a result, simulated discharge and solute breakthrough curves can be compared directly to the actual quantity and water quality measurements taken at the Dickenberg adit.

Hydrogeochemical Conditions
All scenarios simulated in this work included isothermal conditions at 25 °C, with air and liquid as phases subject to transport. The air was composed exclusively of gaseous oxygen as it is the only relevant gas in the oxidative processes. The oxygen behaviour was assumed to follow that of an ideal gas and to be in equilibrium with the aqueous solution throughout the simulation. The chemical composition of the liquid phase corresponds to actual formation, rebound, and precipitation waters taken directly from the field or surrounding areas (Table 2). Their use, however, depends on the simulation stage of the model (see "Simulation Stages" section). Other solutes resulting from dilution and diffusion processes between these water types and atmospheric oxygen were also considered throughout the simulation.
The solid phase of the model was constructed based on detailed mineralogical information (Bedoya-Gonzalez et al. 2021a). The Westfield overburden consists mainly of sandstone layers made of quartz, chert, mudstone lithic fragments, muscovite, and varying quantities of pyrite. Although the rocks lack cement, they contain high amounts of authigenic kaolinite as a matrix. There is no report of carbonate minerals such as calcite or siderite. Based on this description, Table 3 lists the abundance of primary minerals assigned to the matrix continuum along with their chemical formulas and dissociation stoichiometry. Regarding the fractured continuum, the primary mineral fraction was assumed to be zero considering its recent origin and element-free volume.
From the mineralogical composition, the dominance of interactions involving pyrite oxidation, kaolinite dissolution, and dissolved oxygen is clear (see the "Geochemical Processes" section for the details of these reactions). Both pyrite oxidation and kaolinite dissolution are kinetically controlled reactions in the model, while the formation of amorphous iron hydroxides proceeds on an equilibrium basis due to its short formation time compared to the simulation periods. Further transformations of solid iron phases into goethite, ferrihydrite, or schwertmannite are beyond the scope of this study and were not modelled.
Although it has not been identified in the area, gypsum precipitation was introduced as a potential sulfate sink given the high anion concentrations reported at the drainage. Thermodynamic and kinetic rate constants are summarized in Table 3. In the absence of in situ kinetic measurements, representative bibliographic values of the site conditions are used here, making them neither exhaustive nor unique. Finally, minerals' reactive surface area is set at 3140 m 2 /m 3 -medium except for kaolinite and muscovite,  1 3 which is assumed one order of magnitude larger due to their plate habit. This surface area is calculated using a four-sphere cubic-closed packing arrangement of medium sand-size grains, which is a plausible simplification for tight rocks due to its packing efficiency of 74%.

Simulation Stages
The geochemical signature of the Westfield drainage is the result of cumulative natural and mining processes that have modified the system over time. Since mining is carried out in several stages, the model had to address variable initial and boundary conditions. This section describes the fivestage modelling framework used to simulate long-term and seasonal changes of the groundwater chemistry, rocks saturation, and discharge volumes in the study area. The starting time step interval at each stage was established as 1 s and permitted to increase to a maximum of 1000 s. The time step sizes were automatically controlled to maintain a Courant number less than 0.5. This stability criterion generates more accurate transport component results, which are then transferred to the reaction calculations. In addition, the criterion reduces the accumulation of errors in geochemical speciation when the results of one stage are used as initial conditions for the next one. A sketch showing the five-stages simulation can be seen in Fig. 4.

Stage 1: Saturation and Composition of the Formation Water
Full water saturation was assigned to the entire mesh to allow each element to equilibrate by gravity and capillarity with a water table located at −80 m.a.s.l. This depth corresponds to the second exploited coal level down to which the regional water table dropped after the Dickenberg seam was mined (a process that fractured the shallow overburden).
The simulation was carried out for 50 years, equivalent to the time the mine was under an intensive mining regime. Subsequently, the chemical composition of the resulting water volume was assigned equal to the water extracted at the mining faces (i.e. at −80 m.a.s.l.-see Table 2). However, as the residual water is localized in the unsaturated zone, it was equilibrated with an oxygen partial pressure of 0.2 bar before starting the second transport phase. Results of this stage gave formation water volumes in equilibrium with the lowered water table and with compositions that are weakly influenced by meteoric water, as one would expect from a fracture-drained matrix.

Stage 2: Chemistry of Groundwater During Mining Operation
The unsaturated overburden was equilibrated with percolated water for 50 years, the elapsed time between the exploitation of the Dickenberg seam (i.e. the lower boundary of the model) and mine closure. For this purpose, the water table was kept at −80 m.a.s.l., while a fixed injection rate of 7.5 × 10 -6 kg/s/m 2 (≈240 mm/yr) was applied to the first element of each column. This rate corresponds to the average annual recharge of the area, equivalent to 1/3 of the precipitation. As for the injected water composition, it was calculated by applying evaporation to the rainwater data reported at the French-German border (Herrmann et al. 2006) until the remaining water volume matched one third of the original volume (Table 2). Dissolved oxygen was included in this water with a fixed concentration of 10 mg/L. It should be noted that oxygen content normally varies from 0 to 8 mg/L in rainwater depending on a variety of factors, such as atmospheric conditions, temperature, and the presence of other dissolved substances (Chughtai et al. 2014;Kasmin et al. 2016;Komabayasi 1959). The concentration assumed here was slightly higher than the reported upper limit for simple rainwater, but in line with aeration processes in the first meters of unsaturated zone and mixing with freshwater sources at the ground, both processes expected in the area and that increase the dissolved oxygen of the infiltrated water (Hem 1985;Marcy et al. 2022;Pawlicka 2003). Moreover, constant diffusion and advection of gaseous oxygen from the atmospheric boundary was enabled. At the end of the stage, an equilibrium between percolated and formation waters was reached, which includes dilution of the latter. Results also depicted iron, sulfate, and hydrogen concentrations generated by the rock-water interaction during mining.

Stage 3: Groundwater Rebound
In this stage, the lower free drainage boundary was changed to fully water-saturated, with a constant pressure of 1 atmosphere. This shift allowed the simulation of the gravity-capillary interaction between the unsaturated sequence above the adit and the newly imposed water table during the rebound process. It also served to evaluate the potential exchange of pollutants brought by the mine waters. The rising water would tend to dissolve highly soluble efflorescent salts (a solid by-product of sulfide oxidation) accumulated over decades, flushing them out in much higher concentrations than those measured during active mining (Gzyl and Banks 2007;Younger 1997). This process, usually referred to as the first flush, results in environmentally aggressive waters, with very low pH and high metal loads. In this step, the composition of the first flush corresponded to actual measurements recorded in the decommissioning plan of the Ibbenbüren mine ( Table 2). The mean annual water composition of the first-year flush was imposed on the water-saturated lower boundary. As the water came from the saturated zone, an oxygen partial pressure of 1 × 10 -70 bar was assigned to it, converting it into reducing water. Finally, a conservative artificial tracer was added to the water to better quantify its movement through the overburden.

Stage 4: Post-Mining Water Evolution
After rebound, the concentration of contaminants decreased exponentially until they reached levels close to the natural values observed in the surrounding aquifers Younger 2000). Therefore, the fourth stage simulated the geochemical evolution of the water discharged from the unsaturated column after reaching capillary-gravity equilibrium. This was done by injecting the average annual recharge rate for another 35 years directly into the first weathered element and letting it percolate through the columns ( Table 2). The lower boundary was returned to an unsaturated condition enabling free gravity drainage of percolating water into the Dickenberg adit. At the end, the simulated breakthrough curves were compared with the real measurements recorded in the Ibbenbüren decommissioning plan (Klinger et al. 2019) to assess the accuracy of the model.

Stage 5: Seasonal Variation of the Water Composition
Slight fluctuations of the extracted water composition were observed in quarterly chemical measurements of the years

Results
Prior to reactive transport simulations, the flow component of the MINC model was tested independently to spot potential errors in using the parameters and configuration from the dual continuum model. The results showed a large discrepancy with respect to the measured discharges for 2016 and 2017 when using the same calibrated values (Fig. 5-green line). The recession limbs of the discharge curve indicate greater suction and water storage in the matrix continua during winter months rather than discharge through the fractured continuum. Therefore, it was necessary to reduce the α value of the Mualem-van Genuchten function for the fractured continuum from 0.2 to 0.1. This lower value yielded higher water retentions in the fractures, counteracting the large capillary suction generated by the inclusion of a small matrix continuum next to the fracture. Once this change was made, simulations of the MINC model exhibited good agreement of both measured and simulated drainage discharges with the DC model ( Fig. 5-orange dotted line). The latter implies that much of the discretization, parameterization, and recharge values used in the DCM are still valid for simulating the flow component of the multiple continua reactive transport model.

Benchmark Scenario
Stage 1 Figure 6 shows water saturations for one of the model's representative columns after 50 years of dewatering. It is seen that sandstones and shales layers remain highly saturated, with values ranging between 50% and 90%. Conversely, the weathered elements and the fractured continuum display water saturations closer to their residual amounts (Table 1). These results, obtained from the equilibrium of the sequence with a lowered water table, are relevant to avoid potential over-or under-estimations in the reactive and transport processes included in the model. For example, it is not the same to equilibrate pyrite with oxygen when there is an initial air saturation in the sandstones of 66% (as would be the case if the irreducible water content was taken as the initial saturation of the system) than when there was only 30%. The same sensitivity problem would occur with the volume of saline formation water that has to be diluted. As mentioned above, water and gas saturations calculated in this step were used as initial conditions for the second stage.

Stage 2
The conservative behaviour of chloride was used to evaluate the transport of solutes initially present in the matrix throughout mining (e.g. formation water). During the first 10 years, the initial chloride content in the sandstones (3700 mg/L) was actively washed out, generating discharges with concentrations up to 450 mg/L (Fig. 7a). This steep increase was the result of solute exchange from the first two matrix continua closest to the fracture, as illustrated in Fig. 7b (Missing or smaller blue and orange bars for the 5-and 10-year periods in each layer). Beyond the 10th year of simulation, the chloride concentration in the discharge was significantly reduced due to the constant solute dilution in these two continua. However, the reduction slowed down between year 15 and 20 due to a sustained input of the anion from the third continuum. For the last three decades, chloride discharge concentration continued to gradually decrease as the interaction front of the meteoric water with the formation water moved toward and into the fourth matrix continuum (i.e. rock volumes more than 1.5 m from the fracture). Besides horizontal transport, the model further shows how chloride was extracted first in the layers closest to the surface, generating a vertical dilution front with the underlying layers of the same continuum (Fig. 7b). This phenomenon occurred due to the vertical advection of meteoric water into the matrix, which was also facilitated by horizontal exchange with the fractured continuum. At the end of the 50-year period, the simulated chloride discharge of 180 mg/L agreed well with the 200 mg/L recorded in mining reports of the time (Bässler 1970). For this period, chloride dilutions of up to 90% and up to 30% were also observed for the first two continua of the sandstone layers and the fourth continuum of the upper strata, respectively (Fig. 7c).
Stage 2 further returned iron and sulfate discharge concentrations consistent with historical values measured at the Ibbenbüren Westfield. The simulations of pyrite oxidation produced progressive increases of both ions in the discharge until reaching nearly steady-state levels of 1100 mg/L for sulfate and 250 mg/L for iron in the 20th year (Fig. 8a). When plotting the pyrite distribution in the columns, we observed that the mineral dissolved almost homogeneously over time up to 1.5 m next to the fracture continuum. This pattern was generated by the limited advection and diffusion of gaseous and dissolved oxygen into the matrix once it was exchanged from the fracture. Likewise, the rate and amount of oxygen consumed by the oxidative reaction in this interval caused pyrite dissolution to decrease by half in the fourth continuum (Fig. 8c, d). The latter caused the production of both ions to nearly stabilize over the last 30 years as no further amount oxygen could be introduced into the fourth continuum nor could the solute produced there be transported more rapidly back to the fractured continuum.
Oxygen distribution along the continua was, in turn, related to the development of weathering fronts on both sides of the fractures. As the amount of oxygen decreased in the horizontal direction into the matrix due to pyrite oxidation, less oxygen was available for precipitation of iron hydroxides. Figure 8e shows that the first two continua of the matrix (i.e. the first 50 cm next to the fracture) ends with some amount of the solid iron phase. Most of it precipitates in the fracture continuum where the oxygen content is almost twice as high as in the second and third matrix continuum and 10 times higher than in the fourth (Fig. 8d). Under this scheme, total dissolved iron and sulfate concentrations slowly increase in the solutes of the porous medium, while amorphous iron hydroxides persist to accumulate in discrete sections of the fractured continuum. The latter somehow seems to be contrary to what has been reported in drill cores from the area (Bedoya-Gonzalez et al. 2021a).
Finally, this stage illustrated how pyrite oxidation transformed the meteoric water feeding the Westfield into a strongly acidic solution during its seepage path. In fact, the influence of surface water on the pH of the discharge was erased before the fifth year of simulation. Thereafter, the pH achieved a steady state about 3.2 because of the weak acidity buffering capacity of the kaolinite (Fig. 9a). Dissolution of the mineral occurs in the three continua closest to the fractures, with dissolution rates much lower than those recorded for pyrite (Fig. 9b). When kaolinite was removed in an alternative scenario, a slight increase in hydrogen concentration was observed during the first five years of simulation, so low pH values were reached Fig. 6 Water saturation in one of the model columns; for subsequent column graphs, layers G, K, and L correspond to shales more quickly (Fig. 9a). It is also noted that pH continued to decrease slowly during the last years. Despite the latter, the difference in pH value at the end of the 50-years simulation was only 0.1 between the two scenarios. Huertas et al. (1998) reported kaolinite dissolution rates of about 5 × 10 -12 under atmospheric conditions. This low rate results insufficient to affect either the neutralization of groundwater acidity in the area or the generation of excessive amounts of aluminium and silica in mine drainage due to mineral dissolution.

Stage 3
The influence of the mine water rebound in the unsaturated sequence above the adit was simulated by placing a watersaturated Dirichlet boundary condition at the base of the domain. This boundary included a conservative tracer to quantify the advective exchange of substances from the imposed water table. The simulations revealed that the capillary action of the unsaturated sequence sucks the rebounded water during the first 840 days. For the artificial tracer, this represented a total exchange of 50 times the initial imposed amount on the lower boundary per square meter (Fig. 10a). Most of this was exchanged during the first 100 days, with values of about 35% per day. However, water suction progressively decreased as the sequence became saturated. By day 840, the overburden reached hydrostatic equilibrium with the air-entry pressure and the recharged water, stopping the capillary rise. As for transport, the solute moved 10 m toward the surface within the overburden, reaching concentrations one to three orders of magnitude less than what was originally present (Fig. 10b). This thickness would correspond to the capillary fringe where solutes remain above the water table against gravity. The simulation of an alternative scenario with a homogeneous sandstone overburden, however, showed that the thickness of this fringe and the amount of exchanged tracer depend on the rock texture. In the latter scenario, upward migration increased up to 20 m (Fig. 10c), with tracer exchanges about 85 times the initial amount per square meter. Similar capillary patterns to those obtained for the tracer were observed for the sulfate and iron ions. Recall that concentrations of 4500 and 2000 mg/L, respectively, were imposed on the water table boundary to emulate the large amounts measured in the mine drainage during the first flush. In this case, the simulations showed an enrichment of dissolved iron and sulfate for the solutes contained in the two lowermost layers of the model (Fig. 10d). In none of the simulations did the pore water reach saturation levels with respect to gypsum, nor was there an increase in iron hydroxide precipitation that would promote the removal of any of the dissolved ions.

Stage 4
With gravity-capillary equilibrium achieved, water flow was reversed from the column to the water table. At this point, the lower boundary of the model was desaturated to simulate an open mine drainage, avoiding the re-suction of the solute. Figure 11a shows the evolution of the remaining tracer for both the current heterogeneous overburden and a hypothetical rock sequence composed only of sandstones. In both cases, more than 70% of the tracer was released back into the adit during the first 10 years after capillary equilibration. From this point on, the release slowed down as the meteoric water had to travel longer distances to interact and dilute the solute of the inner matrix regions. The process was additionally counteracted by the diffusive movement of solute toward the surface. For example, Fig. 11b shows tracer signals one layer higher than at the end of the capillary-gravity equilibrium for the actual heterogeneous overburden (Fig. 10b). At the end of the 35-year simulation (i.e. year 2018), 20% of the initial amount of absorbed tracer remained in the actual heterogeneous overburden, while only 10% was present in the homogeneous sandstone overburden.
For iron and sulfate, the post rebound stage was also marked by the dilution and release of the absorbed quantities. At the end of the capillary equilibrium (i.e. 1985), sulfate and iron discharge concentrations of the unsaturated sequence were of the order of 2000 and 450 mg/L, respectively. During the following four years, concentrations of both ions decreased by 25% following the drainage of the two matrix continua next to the fracture (Fig. 12a,  b). At that time, concentrations kept decreasing, but only slightly, due to the lower contribution of solute contained in the more distant continua, as was observed for the tracer.
The simulated post-rebound concentrations have similar discharge patterns to those measured, although with lower values. Model underestimations during the first 10 years would arise from not accounting for dissolution and transport of fluorescent salts directly into the adit. Another uncertainty may be due to higher-than-assumed sulfate and iron contents in the rebounded water that was brought into contact with the unsaturated sequence for its capillary equilibrium. In any case, it is observed that the external source was exhausted over time, making the simulated values progressively closer to the measured ones. After 2003, both simulated and measured sulfate signals reached a quasi-steady state close to the pre-rebound discharges, suggesting that the ion concentrations again depend on the water-rock interaction above the drainage (Fig. 12a). Accordingly, the observed difference between the reference levels from this point on could be explained by underestimation of parameters already considered in the model, such as pyrite amount, reactive surface area, or fracture density. Iron discharge behaviour, however, is different from that of sulfate since it did not reach a stable level but further decreased with time (Fig. 12b). A similar decreasing pattern was obtained with the model by including a lower boundary with gaseous oxygen to emulate an aerated adit. This boundary condition leads to the precipitation of iron hydroxides in the sandstone layer immediately above, thus reducing the free iron discharge (Fig. 12c).

Stage 5
The model´s effectiveness in simulating seasonal discharges of sulfate, chloride, and iron was evaluated by comparing quarterly measurements taken at the outflow of the Dickenberg adit. Figure 13a and Fig. 13b show positive agreements of the anion patterns for the years 2016 and 2017, with some discrepancies in the absolute concentration values, as evidenced in the post-rebound stage. Specifically, high recharge events led to lower sulfate discharge concentrations due to dilution and less availability of gaseous oxygen to oxidize pyrite. On the other hand, periods of low infiltration led to higher sulfate concentrations for matrix solutes, which were slowly released to fractures. For chloride, high periods of precipitation generated a slight increase in concentration as there was an increase in the amount of water exchanged between the fracture and matrix continua. This interaction caused more of the solute locked into the matrix to be pushed to the fractured continuum. In dry periods, when the recharge was minimal, the exchange was restricted, keeping chloride in the matrix. Regarding iron, the simulated concentrations agreed with the measured base level, but not with the discharge pattern (Fig. 13c). The measured cation concentrations were shown to be rather independent of the recharge events, which is contrary to the simulated signal, where peaks and valleys were similar to those observed for the anions.

Variations to the Model Setup (Case Scenarios)
This study simulated four additional scenarios to determine the effect of some uncertain parameters on the discrepancy between simulated and measured ion baselines. For this exercise, sulfate ion concentrations were compared, due to their conservative behavior in the absence of gypsum precipitation and bacterial sulfate reduction. Each scenario started from the benchmark setting explained in the methodology section, incorporating the parameter's local variation. All scenarios then went through steps 1 and 2 described in the simulation section to obtain the discharged sulfate concentration for the Dickenberg adit. The details of each case scenario are summarized in Table 4.

Case Scenario 1: Grid with Two Continua (DCM)
This scenario evaluated the dual continuum grid setup employed in Bedoya-Gonzalez et al. (2022). The comparison of the discharged sulfate between the benchmark scenario and the DCM configuration showed strong sensitivity to the number of employed continua (Fig. 14a). When only one continuum was used for the matrix (i.e. the DC model), pyrite dissolution was strongly underestimated as sulfate concentrations increased much slower than when the number of continua was increased to four. The higher anion production of the MINC setup at earlier times was generated by making the continua spacing small enough to resolve horizontal advection and diffusion gradients of gaseous and dissolved oxygen around the fracture-matrix interface. Only towards the end of the 50-year simulation did the DCM show relatively good approximation to the benchmark baseline, when production and discharge of sulfate depended on the large volumes of rock in the innermost part of the matrix. The latter process is evident when looking at similar amounts of average pyrite dissolution of both models (Fig. 14b).

Case Scenario 2: Higher Fracture Density
To the authors' knowledge, no study has dealt with fracture density, so continua with 2 and 3 fractures per block were arbitrarily assessed. In both cases, the simulations showed modest variations for sulfate production, neither exceeding 15% of the Benchmark setup (Fig. 15a). Likewise, congruencies in pyrite dissolution were observed throughout the sequence among the three scenarios (Fig. 15b). The highest values were obtained with the three-fractures/block scenario, where more of the matrix was in contact with the oxygenated water, generating greater pyrite oxidation and sulfate discharges. Conversely, the number of fractures was found to be inversely proportional to the precipitation of iron hydroxide. Higher fracture densities made more oxygen available for pyrite oxidation, generating less for precipitating iron hydroxide. Figure 15c shows that precipitation peaked mostly in the fractured continuum above the shale layers, where less seepage water was absorbed, and thus, more oxygen was available to carry out the precipitation process.

Case Scenario 3: Pyrite Oxidation Kinetics and Reactive Surface Area
Given the lower values reported in the literature (e.g. Nordstrom and Alpers 1997; Williamson and Rimstidt 1994;Xu et al. 2000), this scenario reduced the kinetic rate constant for pyrite oxidation by an order of magnitude. The result was a notable reduction of discharged sulfate compared with the benchmark scenario (Fig. 16a). This indicates that the reaction rate is a sensible parameter in the Westfield primarily  Table 4 Specifications of the different flow, transport and reaction parameters applied to each case scenario Case ID Changed feature from the benchmark setup Observation/reason 1 Grid simplification to 2 continua Employment of the original grid discretization used for modeling the flow component of the Westfield (equivalent to a dual continuum model) 2 Fracture density increases to 2 and 3 fractures per 10 m length block. Kinetic rate for pyrite oxidation is decreased by 1 order of magnitude (log k = 1.0 × 10 -10 mol m −2 s −1 ) To the best of the authors' knowledge, there is no study on fracture density derived from coal mining 3 Reduction of the reactive surface area of the minerals base on a simple cubic packing with eight spheres of medium sand grains ( A s =1600 m 2 /m 3 of medium) Due to lack of in-situ measurements, lower kinetic and surface area values are tested. Kinetic rate are taken from Xu et al. (2000) 4 Variable pyrite volume in sandstones and anthropogenic waste rock deposits (see supplemental Table S Comparison of scenarios with 1 (benchmark), 2, and 3 fracture densities for a sulfate discharge evolution, b pyrite oxidation (i.e. dissolution), and c iron hydroxide precipitation in one of the model's columns because of the constant oxygen supply from the atmospheric boundary. When the latter was turned off, the sulfate concentration was almost equal to that assumed for the recharged water since little pyrite can be dissolved when oxygen is limited to that dissolved in the recharge water.
As there was also some uncertainty about the reactive surface area to be used, the present scenario alternatively reduced the reactive surface area of the minerals to 1600 m 2 /m 3 of material. The latter value was calculated using a simple cubic packing with eight spheres formed of mediumsized sand grains. The simulation results reveal a decrease in the amount of discharged sulfate, although not in a proportional way. While the surface area was reduced by half with respect to the benchmark, the sulfate concentration was only reduced by 25% (Fig. 16a). This disproportion may indicate a co-dependent sensitivity of the parameter, probably linked to the exchange of water and oxygen between the fractured and matrix continua. The latter would be supported by the lower amount of dissolved pyrite at the end of the 50-year simulation (Fig. 16b).

Case Scenario 4: Variable Amounts of Pyrite in Waste Rock Deposits and Sandstones
As mentioned in the setup section, this model included materials from two different origins: (1) the Carboniferous overburden, composed of sandstone and shale layers; and (2) waste rock deposits derived from mining. In this scenario, the pyrite content in these two domains were changed independently. The sole discharge of the waste deposits ( Fig. 17-markers only) showed strong increases in sulfate concentrations when their initial pyrite content was increased. This would indicate that neither the amount of oxygen nor the oxidation kinetic rate was limiting factors in this zone for pyrite dissolution. The same pyrite increases, however, generated minimal variations in sulfate concentrations at the Dickenberg adit, where water from the entire area was collected (Fig. 17-solid lines with markers). The sparse areal representation of the waste deposits in the Westfield produced a difference of less than 10% in sulfate concentrations when comparing the benchmark simulation to the scenario with 30% volumetric pyrite.
The influence of pyrite variations within the Carboniferous sequence was evaluated through three sub-scenarios. In each of them, the sandstones' pyrite content was randomly varied within ± 2 standard deviation from the counts reported in Bedoya-Gonzalez et al. (2021a), while maintaining a mean of 5% for the entire column. A fourth sub-scenario with 10% pyrite content per layer was also evaluated in order to include the most extreme case, with the highest amount Benchmark (Py=5%) Only waste rocks discharge Py waste rocks=20% Only waste rocks discharge Py waste rocks=30% Only waste rocks discharge Fig. 17 Comparison of sulfate discharge for scenarios with varying amounts of pyrite in waste rock deposits. The graph shows concentrations recorded at the base of the waste rock deposits (markers only) and in the Dickenberg adit, where water from the entire area was collected (solid lines with markers) 1 3 of pyrite recorded for the rocks. Details of the amount and distribution of pyrite per layer for each sub-scenario can be found in supplemental Table S-3. Figure 18a shows an overall model insensitivity to heterogeneous pyrite distribution. The opposite occurred when pyrite content was homogeneously increased for the whole area. A doubling of the mineral volume generated a significant increase in discharged sulfate. As with waste rock deposits, the limiting factor for iron and sulfate production was deemed to be pyrite availability within the first two matrix continua rather than the amount of dissolved oxygen. In fact, oxygen surplus causes iron hydroxide to still precipitate in the sub-scenario with 10% pyrite, although its amount decreased greatly compared to the benchmark setup (Fig. 18c). Finally, Fig. 18b, c depict how the heterogeneous distribution of pyrite in subscenario 1 shifts the dissolution and precipitation peaks repeatedly observed in the previous simulations, although maintaining the same base values. This was also true for sub-scenarios 2 and 3, which were not included in these graphs to improve their readability.

Discussion
The construction of a multiple continua model from a dual continuum setup, which in principle is easier to parameterize, proved to effectively simulate solute transport in the unsaturated overburden of the Westfield. Only the α parameter of the fracture continuum had to be reduced to counteract the strong gradient that arises when shortening the mesh node distances with the new matrix elements. The latter caused intense fluid mobilization from the fracture to the matrix, in contrast to reality, as percolated water tends to flow mostly through the fracture network. It is important to note that van Genuchten's α parameter was inversely related to water intake suction for wetting events (Benson et al. 2014). This means that the lower the α value, the greater the water retention within the fractured continuum. Although the capillary parameters for fracture networks are not well understood, the values used here are consistent with the idea of imposing some that reflect the behaviour of a gravel-like medium, as suggested by Bedoya-Gonzalez et al. (2022). Geochemical simulations further revealed the necessity of extending the DC model to a MINC. It was observed that a model with only two continua cannot reproduce the exchange and movement of oxygen within the matrix, greatly reducing pyrite oxidation in the overburden. Similar drawbacks when characterizing the chemical component of a fractured system using the dual continuum approach were reported by Lichtner (2000) and Iraola et al. (2019). Both studies revealed under-representations of kinetic reactions with short equilibrium length scales, as well as numerical instabilities due to the lack of chemical and flow gradients within the matrix. Thus, dividing each matrix block into a finer grid of nested elements seems to overcome the previously reported limitations. Increasing the number of continua results in greater accuracy when calculating diffusive and advective exchange for a system with two high contrast domains. Results from the benchmark scenario in this study certainly provide a realistic geochemical evolution of the system, even accounting for interaction fronts for pyrite oxidation and iron hydroxide precipitation, as described by Bedoya-González et al. (2020) on drill core samples from the Westfield.
Besides grid extension, the model required the use of a five-stage approach to address the limited subsurface information at the beginning of the mining period many decades ago, and the modification of boundary conditions due to anthropogenic activities. Simulations from stages 1 and 2 yielded saturations and water compositions consistent with the development of a water-conductive fracture zone within a moderately permeable overburden. Specifically, dilution fronts of geogenic solutes (dissolved chloride for this study) were formed in both vertical and horizontal directions due to advective and diffusive exchange between the fracture and the matrix elements. This behavior would reevaluate the concept of a single depth-linked salinity gradient for this type of postmining systems. Indeed, application of transient recharge data in stage 5 of the simulation corroborated the interaction and progressive release of solutes in the horizontal component. Typically, intense recharge events in equivalent systems produce major dilution processes that decrease the ion concentration in the discharge (Lasagna et al. 2013;Robins 1998;Whittemore et al. 1989). In the Westfield, chloride concentration showed the opposite behaviour. Higher recharge events caused greater horizontal displacement of meteoric water into the matrix, drawing the locked chloride back into the fractured continuum shortly after. However, over the years, the concentration of chloride in the discharge does tend to decrease due to the depletion of the solute in the matrix.
Results from stages 1 and 2 also showed acidic discharges, characterized by high iron and sulfate concentrations, which are congruent with measurements taken at the Dickenberg adit. This agreement was produced from constant pyrite oxidation in the unsaturated zone and promoted by diffusive oxygen supply from the surface (here represented with an atmospheric boundary condition). The relative thinness of the overburden led to constant flow of oxygen through the vertical fractured continuum and horizontal exchange with the matrix. The latter process was limited to distances of up to 1.5 m from the fractures, mainly due to two factors: (1) matrix permeability restricting dissolved oxygen advection and gas diffusivity and (2) the oxidation rate of pyrite consuming the gas. Therefore, the stable ion discharges achieved after 30 years of mining are mostly derived from pyrite dissolution in this rock volume. It was also seen that the amount of gaseous oxygen in this interval governs sulfate variations when time-dependent recharge was applied in the fifth stage. For instance, intense recharge events generated higher water exchanges from the fracture to the matrix, reducing the amount of oxygen available to dissolve pyrite. In the case of iron, a flat concentration trend that does not react to variable recharge events was observed. This is attributed to a more complex production of hydroxide from oxygen gas at the lower boundary (i.e. at the discharge point to the adit) than the one considered in this study.
Alternatively to the geogenic interactions, the third and fourth simulation stages calculated the influence of the rebound process on the compositions of pore and discharge waters in the unsaturated sequence due to capillary equilibrium. In both stages, the lithological arrangement and water saturation of the overburden (obtained from stages 1 and 2 of the model) were relevant parameters in the quantification of exchanged solute and its vertical transport. For the actual heterogeneous overburden, the suction of a conservative tracer from the rebounded water table was 50 times the initially imposed concentrations per square meter, with a transport up to 15 m upward. In the case of a homogeneous system composed exclusively of sandstones, the values increase to 80-fold and 20 m, respectively. Both estimations represent worst-case scenarios, as the model does not take into account phenomena that would reduce capillary action or contaminant loads, such as interruptions in the continuity of the gas phase, reduction of the porous wetting properties, thermodynamic variations, or stratification of the water column by density (e.g. Boduroglu and Bashir 2022;Costanza-Robinson and Brusseau 2002;Hassanizadeh et al. 2002;Mugova and Wolkersdorfer 2022). Thus, the calculated transport distances correspond to the minimum margin necessary to prevent contaminants from reaching the surface or overlying aquifers. In the Westfield, about 96% of the land surface is above the 15 m distance, while 92% would be above 20 m. After 35 years, capillary forces maintain the thicknesses of these intervals although with concentrations of 10 to 20% of the initially exchanged amount owing to the interaction of the percolated water with the matrix continuum. The latter process, in fact, generated vestigial enrichment of sulfate and iron discharges, which mostly dissipates within 20 years after rebound.
Throughout the five-stage simulation, the model was able to additionally describe precipitation of iron hydroxide with some accuracy. Out of the five continua, the fractured and first 0.5 m of the matrix next to it (i.e. the first two matrix continua) presented saturated and oxygenated conditions for their precipitation. This pattern was consistent with the weathering fronts described by Bedoya-González et al. (2020) for the Westfield sandstone layers, with the exception of fracturefilling hydroxides. The latter discrepancy is understandable considering that mineral precipitation in real open fractures is a complex process that depends on the right amount, size, and arrangement of heterogeneities within the medium to seed the particles (Jones and Detwiler 2016). It has also been found that transport limitation resulting from slow flow velocity, deadend fractures, or isolated structures favours the precipitation of solid phases (Deng and Spycher 2019;Menefee et al. 2017;Singurindy and Berkowitz 2005). Then, the extension and high permeability of mining-derived fractures in the Westfield would be an impediment to solid-phase precipitation. The amount of iron hydroxide that simulations revealed to precipitate within them may have actually been transported as suspended particles into the Dickenberg adit.
Finally, results from four alternative scenarios complemented the findings obtained with the benchmark setting. Although all four gave similar pyrite oxidation fronts relative to the fractured continuum, iron and sulfate discharges did vary according to the modified criterion. Based on actual Westfield features, equal ion discharges to those recorded in recent years (i.e. 1600 mg/L sulfate and 200 mg/L iron) could be reached by increasing pyrite oxidation rates by less than an order of magnitude or a density of 3 fractures per 10 m. For the first case, the increase could be derived by the action of iron-oxidizing bacteria, as has been widely discussed in the literature (e.g. see Johnson 2003;Johnson and Hallberg 2005;Nordstrom 1982Nordstrom , 2011. Fracture density, on the other hand, is a more complicated parameter to evaluate since, to the authors' knowledge, there is no study that makes an approximation of it. However, Bedoya-Gonzalez et al. (2022) showed that this fracture density does not notably affect the flow behaviour of the Westfield, so both chemical and transport simulated components would be consistent with what was measured. Pyrite reactive surface area was also a sensitive parameter, but the value used for the benchmark setting was already high and consistent with the lithologic and hydraulic conditions. Moreover, the use of geometric relationships seems to be a good approximation for regional models due to uncertainties that arise when scaling measurements taken under controlled conditions of pH, temperature, and grain size. Lastly, the restricted location of the waste rock deposits makes the volume of pyrite within them insensitive. The same is true for the heterogeneous spatial distribution of the mineral along the layers as long as the average content in the overburden remains at 5%. However, varying the average pyrite content of the entire sequence does have a marked effect on the ion load discharge.

Conclusions
We formulated and applied a multiple interactive continua (MINC) model to characterize advective and diffusive reactive transport processes in fractured post-mining coal zones. The employed approach represents a novel attempt to simulate water-rock interaction compared to commonly employed lumped or equivalent continuum models, as it simultaneously recognizes fracture and porous flow regimes. Its effectiveness was evaluated here by simulating the origin, generation, and transport of AMD within the shallow and fractured overburden of the Ibbenbüren Westfield. The model was constructed by extending a dual continuum model, previously used to calculate the flow component of the zone, to a five-continua model (one for the fracture and four for the matrix) that accounts for reaction fronts and diffusion through the matrix blocks. Additionally, dissolution, precipitation, and dilution processes linked to iron, sulfate, and chloride-bearing minerals and solutions were included using a five-stage simulation process. The latter was designed to address variations in geochemical and hydraulic conditions derived from the mining history of the area.
Simulation results produced discharges with ion concentrations close to those measured in the Westfield. Horizontal interaction between the formation water contained in the matrix continuum with the water percolating through the fractures would be responsible for the high chloride concentrations measured in the mine adit. This process also results in dilution fronts in both vertical and horizontal directions for solutes initially present in the matrix, generating quasistationary discharges over time. Complementary flow of gaseous oxygen from the atmosphere through the fracture continuum proved to be highly relevant in the study area. Subsequent horizontal diffusive exchange of oxygen with the matrix continua closer to the fractures generates constant pyrite oxidation that releases acidity, sulfate, and iron loads to the water. This oxidative process alone would produce the long-term and seasonal concentrations measured in the Westfield, if small modifications in fracture density, oxidation kinetics, or mineral surface area are included, none of which affect the flow component of the model. The exception to the good fit between measured and simulated concentrations was the elevated iron and sulfate concentrations measured during the first 10 years after rebound. The difference was evaluated here by simulating the capillary suction of the rebounded mine water from the unsaturated sequence along with its subsequent release. Although the measured values did not match the simulated values, the coincidence in the decrease of the concentrations of both curves corroborate the existence of a capillary equilibrium between the unsaturated sequence and the rebounded mine water. The latter was used to establish safety levels in the unsaturated overburden by determining maximum concentrations and vertical transport distances of the rebound water during the creation of a capillary fringe. It was also used to predict the concentrations of long-lasting solutes in the capillary fringe decades after the rebound process.
Acknowledgements The authors gratefully acknowledge the assistance of Professor Dr. Christoph von Hagke, Dr. Moritz Liebl, and Nicolas Villamizar M.Sc. from the Department of Environment & Biodiversity at the University of Salzburg for their advice with data management and manuscript construction. We also thank the team at RAG Anthrazit Ibbenbüren for facilitating discharge measurements and providing mine details. Our appreciation also goes to the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen (LANUV) for supplying the study area's precipitation data. Finally, we appreciate the comments and suggestions of the two anonymous reviewers and this journal's editorial staff, who helped us to improve the manuscript's quality.
Availability of data and material Input and output files generated for and from TOUGHREACT are available from the corresponding author on reasonable request. Precipitation data were obtained and used after written communication with the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen (LANUV). The corresponding author is welcome to share them with permission from LANUV and upon reasonable request. All other datasets used and/or analyzed during the current study are included in this published article and its supplementary files.
Funding Open access funding provided by Paris Lodron University of Salzburg. This work was financially supported by Forum Bergbau und Wasser (FBW).
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/.