To study the location, density and the pattern of fractures in the top of the Cretaceous, an ERT survey was conducted. ERT is used to determine the electrical resistivity of geological layers, by achieving 2D or 3D inverted sections (Loke et al., 2013). ERT is suitable for many hydrogeological applications (Binley et al., 2015) such as locating cavities and fractures or estimating layer thicknesses and detecting interfaces between layers (e.g., Guérin, 2005; Al-Tarazi et al., 2006; Guérin et al., 2009; Robert et al., 2011; Van Hoorde et al., 2017). It has also been used together with tracer experiments to identify preferential flow paths and solute transport processes along fractures (e.g., Robert et al., 2012) or image infiltrating processes in the unsaturated zone (Lesparre et al., 2017; Claes et al., 2019; Blazevic et al., 2020).

In this study, we assume that the fractures are creating a secondary porosity in the rock matrix, therefore increasing the water content in the saturated zone. The resistivity of the chalk is much higher than that of water (e.g., Reninger et al., 2014; Katika et al., 2018). In fractured zones, more water is present. According to petrophysical laws (Revil et al., 2017), this is accompanied by a decrease in electrical resistivity. We also expect a contrast in resistivity between the consolidated chalk deposits and the overlying unconsolidated clastic (sand, silt, and clay) sediments.

A total of six profiles were arranged in the studied catchment with different lengths and in different directions traversing both the valleys and the hillsides (Figure 2) to accurately image the local pattern of fractures. The first four profiles were initially laid out to verify if fractured zones were limited to the valleys as originally assumed. Profiles 5 and 6 were later collected to assess the continuity and orientation of the previously identified fracture zones. Four profiles (profiles 1, 3, 5, and 6) were set out around the four pumping wells in the center of the study area, while profiles 2 and 4 are located further away on hill slope and crests.

The ERT data were collected using an ABEM terrameter LS with 5 m electrode spacing. A multi-gradient electrode configuration was adopted (Dahlin and Zhou, 2006), given its relatively good resolution to detect lateral variations and its high signal-to-noise ratio. Each initial profile comprised 64 electrodes and 1,392 measured resistances. However, when possible, the profiles were extended using the roll-along technique.

It is crucial to accurately quantify the noise to prevent misinterpretation of ERT images (LaBrecque et al., 1996). Overestimation of the noise can lead to gross smoothing of structures, while an underestimation of noise can create artificial structures (Slater et al., 2000). Noise can arise from several factors like poor electrode contact, random errors attributed to the measurement device and errors arising from external factors (background noise) (Slater et al., 2000). To estimate the error on the ERT measurements, reciprocal measurements were collected for a subset of the data set. Overall, an increase of the reciprocal error with raising resistance was observed. Profiles located in the valley showed a very low relative reciprocal error, lower than 1%, while the profiles on the hill slopes and crests showed larger reciprocal errors in the range of 2–4%. For the latter, the substrate was hard after a dry period, increasing the electrode contact resistance, what likely explains the larger reciprocal error. Overall, the quality of measurements is sufficient and can be used for the qualitative interpretation of fractured areas. The error levels are also consistent (within the same range) with the final error of the inversion process (see below).

The data processing and inversion was accomplished with the software Res2dinv (e.g., Loke, 2006). We used a blocky model constraint (L1 norm) as the expected resistivity distribution resulting from fractures is expected to yield sharp contrasts in resistivity, with lower resistivities concentrated in fracture zones. A L1 norm is also used on the data to minimize the effect of outliers. Before inversion, the apparent resistivity pseudo-section was plotted to remove obvious outliers. After a first inversion, data points with a root-mean-square error larger than 100% were also removed from the data set, to avoid potential artifacts of inversion.

The sensitivity was used to assess the quality of the inversion. An example is given for profile 1 in Figure 3. It shows the classical decrease of sensitivity with depth and on the side of the section, as well as a lower sensitivity in resistive blocks as electrical current flows preferentially in conductive zones. Additionally, in order to check the quality of the obtained images and to derive to which depth the ERT results are reliable, we determine the depth of investigation index (DOI) using two inversion with reference model equal to 0.1 and 10 times the average apparent resistivity, as suggested by Oldenburg and Li (1999). The calculated DOI index is a value between 0 and 1. Indexes approaching zero mean that both inversions produce the same resistivity value. DOI indexes close to one are in contrast not sensitive to surface data anymore but depend mostly on the reference model. We use a threshold of 0.1 for the interpretation of the resistivity images (e.g., Caterina et al., 2013). An example of the DOI is shown for profile 1 (Figure 3C). It shows that the whole section has DOI values lower than 0.1, except for a zone at the border of the section, meaning that we can interpret the variations of resistivity with confidence down to depth around 80 m. Slightly larger DOI are observed in some resistive blocks and at the interface between the chalk aquifer and the overlaying unconsolidated deposits, but they do not prevent the interpretation of the corresponding interfaces.

Several layers can be identified on the ERT profiles (Figures 4A–F). The interpreted top of the Cretaceous is marked as a white dashed line at the location of the contrast between the resistive chalk and the more conductive overburden. In profiles 2, 3, and 4, the layer above is identified as the Formation of Heers. It has a lower conductivity due to a high clay content in the Member of Orp and the high water content in the relatively well-permeable crumbly marl of the Member of Gelinden (Vandenberghe et al., 2004; De Bast et al., 2013; Deckers et al., 2014). The top of this formation is indicated with a white dotted line. The upper layer is interpreted as the Quaternary. These sediments generally have a higher resistivity than the Formation of Heers but the contrast can be very low which complicates separating them like in profiles 1, 5, and 6. ERT is thus able to roughly identify the different geological layers in the upper part of the model. Those layers are laterally relatively homogeneous, so that that the recharge pattern can be expected to be mostly dominated by the topography of the area and the soil occupation and not by preferential infiltration path.

Interpreted fractured zones are indicated as full red lines above the profile and correspond to lateral decreases of resistivity in the chalk formations. Smaller contrasts in resistivity are indicated with a dashed red line above the profile, and could correspond to fracture zones with a smaller fracture density or water content. It can be concluded that fractured zones in the top of the Cretaceous are not only found in the valleys (profiles 1, 3, 5, and 6) but also in the hillier parts of the study area (profiles 2 and 4). A second important observation is that profile 3, although located in the valley, does not contain clear fractured regions in the chalk, but only slight contrasts in resistivity. A final important observation is that the fractured zones in profile 1, 5, and 6, seem narrower than the ones in profiles 2 and 4. These observations seem to suggest that the identified fractured zones are subparallel to profile 3 and almost perpendicular to profiles 1, 5 and 6, while they intersect profiles 2 and 4 with an angle.

Most of the available boreholes are limited to the center of the study zone (Figures 1, 2) and it was not possible to systematically measure the profiles in the vicinity of the boreholes. In addition, the chalk aquifer is not outcropping in the study area. For these reasons, it is impossible to validate the identified fractures with ground-truth data. Nevertheless, ERT results are coherent with the expected distribution of geological formations, and the DOI and sensitivity indicate that the identified low resistivity zones in the chalk are in the resolved part of the tomograms. We are thus confident these anomalies correspond to fractured zones and/or zones with preferential dissolution in the chalk. The fracture zones do not seem to generate heterogeneities within the upper part of the model. Only one well is co-located with profile 1 (distance of 275 m along the profile, see Figure 2). No specific indication of fractures was found during drilling, but the well seems rather located in an area where the top of the chalk is globally less resistive, maybe indicating more alteration and dissolution in the top of the chalk aquifer. For tapping in the higher hydraulic conductivity zone, the well should have been placed 50 to 100 m further to the North East.

To better grasp the spatial distribution of the fractures, they are set out on a map in Figure 5. Despite the relatively large uncertainty related to the distance between the parallel profiles, it is possible to propose a tentative connection between some of the fractured zones using the profiles in the center of the study area. The fractures are interpreted to be almost parallel to profile 3, oriented along the northwest-southeast direction. This orientation is also parallel to the relatively nearby Roer Graben (Deckers et al., 2018) and its numerous faults in the western part of Germany, southern part of the Netherlands and eastern part of Belgium. This orientation has been previously identified as the main orientation of fractures (Figure 5B) in the Cretaceous formations (DOV, 2019). The extension of the fracture zones beyond the area investigated by ERT is unknown, but the global analysis of main fractures in the Cretaceous formations tend to show that fractured zones crossing the whole study area can be expected. This assumption will be further investigated with the numerical models.

Based on the available geological, hydrogeological and geophysical data, as well as the objectives of the study to predict the average drawdown resulting from pumping, it was chosen to only model the groundwater part of the hydrological system. The aquifer from which groundwater is extracted is semi-confined so that the drawdown will be dominated by processes taking place in the chalk aquifer. In addition, the upper part of the model is relatively homogeneous, and the recharge of the aquifer from the surface is occurring through a series of low permeability layers, so that for the purpose of the study considering an average recharge rate is sufficient. In addition, the available data are mostly limited to the center of the catchment in the exploited aquifer, so that calibrating a fully integrated hydrological model is illusory and not necessary. Nevertheless, the interaction with surface water will be integrated through the boundary conditions.

The water extraction site of Voort is then modeled using Modflow-2005 (McDonald and Harbough, 1988; Harbaugh, 2005). In the horizontal plane, the grid is constructed with an extent of 15,000 by 15,000 m based on the estimated influence radius of the extraction wells and their distance to natural hydrogeological boundaries. The 4 pumping wells at the main extraction site are located close to each other in the middle of the study area (Figure 2). These four pumping wells together have a permit for a pumping rate of 2,190,000 m³/year (De Watergroep, 2014). Another smaller water extraction site with two pumping wells is present in the North of the model area. This extraction site has an average pumping rate of about 650 000 m³/year over the period 2001–2017. The size of the cells is the smallest around the pumping wells with cell size of 25 × 25 m, the cell size increases progressively toward the sides of the model. Each layer of the model contains 260 rows and columns. Vertically, the model is discretized into nine layers based on the hydrogeological units presented in section Geology and hydrogeology of the study site.

A steady-state approach with varying stress-periods is chosen to simulate the evolution of the hydraulic head in the study-area in the period 2001–2017. Such a simplified approach is necessary due to the lack of data to estimate transient recharge patterns and boundary conditions related to seasonal variations. It is also in accordance with the objective of the study to estimate average drawdowns in the catchment. Eighteen stress periods are used, with each stress period corresponding to the yearly average hydraulic head and extraction rate in the period 2001–2017. The first stress period is defined before pumping occurs, to model the groundwater flow in natural conditions without external influences. This time-step is then followed by 17 other stress periods each representing 1 year with constant pumping at an average rate.

The limits of the model were set beyond the estimated radius of influence of the extraction site and chosen to correspond to natural boundaries whenever possible. This allows to limit the influence of the boundary conditions on the model results. Model cells outside of these limits were set as inactive. For the upper layers of the model, the boundary conditions were set up as follows (Figure 6).

In the South of the study area, the hill ridge serves as the southern boundary of the model. It is considered as a groundwater divide with the adjacent catchment. It can thus be characterized by a zero-flux crossing the boundary (Neumann boundary condition). The western boundary follows a spur of the southern hill ridge and is thus also assigned a zero-flux boundary. From this ridge, the boundary continues perpendicularly toward the Melsterbeek (Figure 6). Along this segment, the water flow direction will be mostly directly toward the river, with limited amount of water crossing the boundary and a zero-flux boundary type is therefore assigned. The Melsterbeek and the Bergbeek are two streams forming the rest of the western boundary. Here a specified hydraulic head is assigned, 1 m below the topography, assuming that the aquifer is in equilibrium with the river. Although the water level can change slightly from year to year, these rivers are located far from the extraction area and these variations are assumed negligible. For the northern and eastern boundary, a similar approach following streams and rivers (Bergbeek, Spaasbeek, Herk River, Kleine Herk, and Mombeek) joined by perpendicular segments is followed, assigning, respectively, fixed hydraulic heads and zero-flux boundary conditions (Figure 6). The hill ridge between the Herk and the Mombeek was considered too close to the extraction site to be used as a boundary. From the Source of the Mombeek a spur of the hill ridge is followed until the southern boundary is reached. The boundary conditions remain the same for all stress periods used in the model.

For the Cretaceous layers forming the semi-confined aquifer, there is a regional groundwater flow component from the South-East to North-West, so that other boundary conditions must be implemented. The regional flow is extracted from the primary monitoring network of the Flemish Environmental Agency (VMM) and specified heads are used along the north-western and south-eastern boundaries. The north-eastern and south-western boundaries are parallel to the regional flow and are thus set as zero-flux boundaries.

The model area is located in the catchment of the Herk river between the Mombeek and the Melsterbeek. Several smaller streams also cut into the landscape. The water extraction site is located around the point where the Molenbeek flows into the Herkebeek. The latter flows into the Herk (Figure 6), which in turn flows into the Demer north of the study area. All rivers and streams in the study area flow toward the lower lying areas in the North (Figure 6). Based on observations, the streams within the study area are relatively shallow and the rivers are therefore implemented as head-dependent flux boundaries (river package). This type of boundary allows an exchange between the aquifer and the river, dependent on the head calculated in the aquifer. The magnitude of this flow is regulated by the conductance factor, which describes the resistance of water flow between the river and the adjacent cells. This factor is in function of the thickness and hydraulic conductivity of the riverbed and the width and length of the river in the cell (Harbaugh, 2005).

The different pumping wells were implemented using the well-package. The pumping rates at the water extraction site in Voort vary between the different wells and between different stress periods. Wells 4017-001 overall has the highest pumping rate with an average of 915 m³/day. This is followed by pumping wells 4017-002 and 4017-004 which pump at an average rate of 692 and 669 m³/day, respectively. Lastly, pumping well 4017-003 has the lowest pumping rate with an average of 433 m³/day. North of Voort, in Wellen, another extraction site is present with two pumping wells, which are never active simultaneously. The average pumping rate of this water extraction site is 1,773 m³/day.

Estimating the recharge is difficult because of many local variations. For the Cretaceous aquifer itself, the recharge zone is located outside the study area to the South where the corresponding formations are outcropping, and an incoming flux is determined by the boundary conditions. For the upper layers, the recharge is coming from rainfall. The total estimated yearly recharge is 198 mm/year (Driesen, 1980), calculated as the difference between the precipitation and the annual evaporation. The latter was calculated using the Coutagne formula (Coutagne, 1954). The recharge is homogeneously distributed in the model and kept constant for the different stress-period. There is no evidence that the presence of the fracture zones in the Cretaceous aquifer have an impact on the recharge patterns, likely because the recharge from precipitation is dominated by the topography and the surface deposits, which are relatively homogeneous. Indeed, part of the recharge is immediately drained by small streams and rivers, as well as through small drains and ditches that are present in the study area. The latter are impossible to model individually. A drainage level was therefore implemented in the model. A linear relationship was derived between the elevation and the depth of the water level in the Quaternary aquifers using shallow piezometers available in the study area and was set as the drainage level. These observation wells are available from the monitoring done by De Watergroep around the pumping wells and from the monitoring network from VMM, from which freely accessible data is available. In consequence, the effective recharge to the model is the difference between the recharge from rainfall and the water drained from the aquifer. The latter is governed by the drainage base and the conductance factor set to 100 m²/day. This approach allows us to simplify the process taking place in the unsaturated zone and at the surface and avoid the complexity of a fully integrated model.

Before calibration, the horizontal and vertical hydraulic conductivities for each of the nine layers of the groundwater were estimated (Table 1). The most sensitive parameters are those from the Cretaceous from which water is extracted (layers 8 and 9). Although data from the Flemish database show local variations in the hydraulic parameters for the different layers, all the layers were considered as homogeneous for the calibration, with initial values taken from the literature (e.g., Calembert, 1956; Maréchal, 1993; Wouters and Vandenberghe, 1994; Diez, 1999; Feyen et al., 1999; Corluy et al., 2004; Lebbe and Vandenbohede, 2004; Vandersteen et al., 2014). A constant ratio of 10 between the horizontal and vertical hydraulic conductivity is kept for all further simulations.

For layer 1 to 7, the initial head is set equal to the drainage level. A first simulation is run for the first stress period (no pumping) with this initial head. The computed hydraulic heads are then used as initial head values for subsequent simulations. This allows to already integrate the effect of rivers in the initial head distribution.

In the Cretaceous (layers 8 and 9), the initial heads are based on an existing piezometric map of the Cretaceous aquifer generated by interpolating the available measured hydraulic heads in the Formations of Maastricht and Gulpen using kriging (VMM, 2008). This distribution is in accordance with the boundary conditions.

Several scenarios (Figure 7) are tested to implement the fractured zones detected by ERT in the groundwater model. In the original scenario A, only the main fractured zones observed in the center of the model are implemented. In scenario B, also the fractured zones detected on the crests and hill slopes are integrated. The other scenarios are assuming that the detected fractured zones are present in the whole model domain, even if ERT data were only available in the center. Scenarios C is similar to scenario A but with the fracture zone extending to the edges of the model. Scenario D is similar to scenario C but with another orientation of the fractures. Scenario E is similar to scenario C but with a unique larger fractured zone instead of 3 smaller ones. Scenario F uses an even larger fracture zones. Finally, scenario G has 3 relatively large fractured zones across the whole model corresponding to the groups of fractures detected by ERT. In the fracture zones, the hydraulic conductivity of the cells of layer 8 (Maastricht Formations) is increased in order to accommodate the fact that water is flowing more easily in the fracture zones.

Although adding anisotropy to characterize flow in faults and fractures is a common approach (Bear, 1993; Medici et al., 2021), the hydraulic conductivity of the fracture zones is modeled as isotropic. This approach is justified because fractures are not modeled individually, but globally as wider zones, and because the contrast of hydraulic conductivity between fractured zones and the background combined with the geometry of the fractured zones is sufficient to simulate preferential flow path along the direction of the fracture.

First, a calibration of the homogeneous model was performed to estimate the hydraulic conductivity of the different layers. The best model fitting was obtained with the parameters of Table 2. In Table 3, the results of different calibration attempts are shown in terms of average misfit (difference between observed and calculated hydraulic heads throughout the time series) for selected monitoring wells located both in the Cretaceous and in the upper layers as well as the maximum drawdown observed in the pumping wells. The objectives being to both reproduce correctly the observed distribution of hydraulic head through time and to correctly estimate the cone of depression resulting from pumping activities.

Note that a perfect fit is not expected as the simulations calculate yearly average water levels using average extraction rates while the data set is made of monthly measurements capturing seasonal variations and variations in the extraction rate. In addition, the size of the model is such that the grid size is not sufficient to accurately simulate individual drawdowns as some strong variations are sometimes observed at short distances. Moreover, pumping wells are affected by clogging, so that the measured hydraulic heads are lower than the actual levels in the aquifer.

The results of calibration of the homogeneous case (Table 3, Figure 8) show that it is very difficult to obtain at the same time an acceptable calibration in all observation wells and a realistic prediction of the maximum drawdown at the pumping wells that should be below 6 m. Although an underestimation of the drawdown at the pumping wells is expected with the model given the averaging effect of the grid size and the clogging of the well, the homogeneous model tends to overestimate this drawdown, clearly indicating an underestimation of the hydraulic head. Some trends observed in the piezometers are also not captured by the model (Figure 8), probably because they result from factors external to the extraction site and are not accounted for in the model (change in the boundary conditions or the recharge rate).

The best result is obtained with a horizontal hydraulic conductivity of 5 m/day for the formation of Maastricht, but the model is highly sensitive to variations of this model parameter. An increase in the hydraulic conductivity of the layer allows to better match the drawdown in the wells, but at the expense of the fit in the monitoring wells, making it impossible to fulfill the two above-mentioned objectives with a homogeneous chalk layer.

Following this reasoning, it seemed that implementing fractured zones with a higher hydraulic conductivity should help to better fit the observed data. The different scenarios involving fracture zones were implemented in the model (Table 3). The hydraulic conductivity of the layer was kept as the one from calibration of the homogeneous case (5 m/day) while the fractured zones were assigned a larger hydraulic conductivity. Globally, adding the fractured zones allowed to better balance between the two objectives of fitting the maximum drawdown in the pumping wells while properly estimating the water level in the monitoring wells. Adding the fractures also improved the fit for wells located in the Member of Orp (layer 7), while the Quaternary layer (layer 1) was not really affected (Figure 8).

From scenarios A (Table 3), it is clear that increasing the hydraulic conductivity of the implemented fractured zones allowed to better fit the drawdown in the well. While the fit in the monitoring wells is degrading, this is not as bad as in the homogeneous case, indicating that including fractures in the model is an efficient way to obtain a better fit. The other proposed scenarios for the inclusion of fractures have only a limited impact on the model misfit. This is likely related to the absence of monitoring wells outside of the center of the model where the main fractures zones where detected using ERT. Nevertheless, since fractured zones are likely not limited to the center of the model, we found interesting to study those alternatives. Globally, they allow to further improve the maximum estimation of the drawdown in the pumping wells, without adding a negative effects on the monitoring wells. A good compromise is found for scenario E3, with a unique, wide fractured zone with a hydraulic conductivity 3 times larger than the homogeneous layer (Figure 8), although it remains difficult to select the best solution given the two conflicting objectives. The quality of calibration is considered as sufficient to investigate the impact of new exploitation scenarios on the drawdown in the catchment. The model of scenarios E3 is kept as a basis for predictive purposes.

In the management strategy of De Watergroep, the extraction rate in the catchment should be increased from 2,600 m³/day to about 6,000 m³/day. Two scenarios are envisaged: (1) increasing the pumping rate of the current wells, (2) drilling a new well in the catchment. For the second choice, the options are limited. To limit the investment costs, the only viable solution is to drill a well next to an existing water reservoir owned by de Watergroep (Figure 9). This reservoir is located on the crest along ERT profile 2. On this profile, a fractured zone was detected, with very low resistivity, and a second possible fractured zone was identified with a less significant contrast in resistivity. Since the water reservoir is located on top of the latter, this option was retained in the simulation scenarios. However, given the limited contrast in resistivity, we expect the hydraulic conductivity to be only slightly larger in the fracture zone.

The objective of the prediction is to estimate the hydraulic head in the Cretaceous layer and the maximum drawdown in the Quaternary. The former will determine the new protection zone around the pumping wells, while the latter is significant for geotechnical considerations. Indeed, the Quaternary layer contains peat, and an excessive drawdown could yield some subsidence and differential settlement in the area. In particular, we analyse the impacts of the conceptualization choices on the model prediction. In other words, we assume that a model including the fracture zones detected by ERT is more appropriate, and estimate how an error in the conceptualization propagates in the model prediction. This is justified as the calibration of the model with a homogeneous chalk aquifer yielded an overestimation of the drawdown and the hydraulic gradient, and therefore an overestimation of the settlement in the Quaternary layer.

Four scenarios are tested (Table 4). Scenario 1 considers the current pumping rate. Scenario 2 considers an increase of the pumping rate using the current wells. Scenario 3 considers a new pumping well-ignoring the fractured zones detected by ERT, with the total pumping rate equally distributed over the 5 wells. Finally, scenarios 4a to 4d considers that the new pumping well is located in a fractured zone, as identified by the ERT results, with increasing hydraulic conductivity. The results are summarized in Figure 10 and Table 4.

With the current pumping rate (Scenario 1), the simulated maximum drawdown in the Cretaceous layers is about 6.2 m, and is 0.2 m in the Quaternary layer. When increasing the pumping rate (Scenario 2), those maximum drawdown increase to 14.19 m and 0.47 m, respectively. Allocating 20% of the extraction rate to a new well located on the crest (Scenario 3), reduces the maximum drawdown in the Cretaceous layer to 13.63 and 0.44 m in the Quaternary layer and significantly reduces the area impacted by large drawdowns (Figure 10). The maximum drawdown remains of course located under the current well-battery, and reducing further the observed drawdown would require to allocate a higher extraction rate to the new wells. In consequence, the presence or absence of a fractured zone only has a limited impact on the maximum drawdown in both studied layers in the new pumping well. However, the drawdown in the new pumping well is significantly reduced from 8 to 6.5 m. Therefore, a better repartition of the extraction rate within the fractured zones detected by ERT would be a reasonable management scenario for this catchment in order to limit the maximum drawdowns observed in the aquifer and overlying layers.