Application of a novel cascade-routing and reinfiltration concept with a Voronoi unstructured grid in MODFLOW 6, for an assessment of surface-water/groundwater interactions in a hard-rock catchment (Sardon, Spain)

Integrated hydrological modelling (IHM) can reliably characterize surface-water/groundwater interactions in complex hydrological systems such as hard-rock systems (HRS), located in water-limited environments (WLE). Such HRS-WLE conditions are represented by Sardon catchment (~80 km2) in Spain, where the MODFLOW 6 modelling environment was tested, applying the following improvements as compared to previous works in that catchment: a new conceptual model, driving forces redefined based on remote sensing data, an unstructured Voronoi grid, and, most importantly, a novel cascade-routing and reinfiltration (CRR) concept. In the standard MODFLOW 6, rejected infiltration and groundwater exfiltration have always been considered as sinks (evaporation). However, in reality, that water can not only evaporate but also reinfiltrate back to the subsurface or move as runoff towards drainage water bodies. The CRR improves surface–unsaturated-zone interactions and also surface-water/groundwater interactions. The standard and new capacities of MODFLOW 6 are presented in the transient model of the Sardon catchment, calibrated using 7 years of daily groundwater heads and streamflows. The results showed: the large spatio-temporal variability of the groundwater fluxes, the substantial role of groundwater exfiltration, the low catchment storage, the fast reaction of the water table and streams to rainfall, and the mosaic character of the net recharge. These characteristics are typical for HRS-WLEs with a shallow water table. MODFLOW 6 has many improvements compared to previous MODFLOW versions, so with the proposed CRR concept (still can be improved), the single-environment MODFLOW 6 has modelling capacity comparable with multienvironment IHMs, while being more flexible and more efficient.


Introduction
Aquifers of hard-rock systems (HRSs) have not been given as much attention as other aquifer types, likely due to their low productivity and difficulties in water-well drilling (Singhal and Gupta 2010). However, in many countries, hard-rock groundwater resources, even though small, represent the only freshwater supply for sustaining life. Therefore, there is a continuous need for improving the assessment of groundwater resources in HRSs, to more efficiently extract groundwater from these aquifers, particularly when surfacewater resources are not available, and to improve management of groundwater resources.
The HRSs are characterized by generally low groundwater storage and spatially variable hydraulic conductivity due to their weathered fractured-rock composition and large heterogeneity. In HRSs, groundwater flow occurs mainly due to the secondary porosity, i.e. faults and fractures, and is largely dependent on their density and connectivity. The low aquifer storage and typically shallow impermeable rock basement, so shallow drainage base, imply that the HRSs are characterized by a shallow groundwater table, dense drainage networks and related short groundwater residence time associated with short flow paths (Hassan et al. 2014). The challenge in dealing with the HRSs is their complexity of surface-water/groundwater interactions and the related difficulties in simulation.
Many HRSs are located in water-limited environments (WLEs), described by the aridity index AI = P/PET ≤ 0.65, where P is annual precipitation and PET is annual potential evapotranspiration (Parsons and Abrahams 2009). The WLEs are characterized by: (1) high spatial and temporal variability of precipitation, with typical erratic but intense showers; (2) large spatial and temporal land-cover changes (type and pattern of vegetation); and (3) vulnerability to desertification, groundwater depletion, salinization, soil erosion and nutrients limitation (Newman et al. 2006). In WLE-HRSs, intensive rainfall showers imply intense recharge and abrupt water-table rise, often resulting in groundwater exfiltration to the land surface. As such, the climatic conditions of WLEs enhance the complexity of surface-water/groundwater interactions in HRSs.
Vegetation cover in WLEs has an important role in the dynamics of surface-water/groundwater interactions. The typical vegetation of HRSs in WLEs is the savannah type of woodland mixed with grassland, i.e. sparse occurrence of patchy woody vegetation, evergreen or deciduous, and short period of growing grass. The sparseness of vegetation and the frequent high-intensity storms lead to the concentration of overland flow and the formation of channel networks (Newman et al. 2006). The vegetation-related processes, i.e. canopy interception and transpiration, limit the replenishment of the water resources, and thus influence surfacewater/groundwater interaction. Moreover, some specific tree species and even shrubs, called phreatophytes, can influence surface-water/groundwater interactions through their tap-root groundwater uptake and its eventual redistribution by shallow lateral root system (a process known as hydraulic redistribution), while different species influence these interactions in different ways (Lubczynski 2009).
Several conceptual models were developed for describing the groundwater flow in HRSs, such as the parallel plate model, double porosity model, discrete fracture network model, and equivalent porous medium model (EPM). The EPM replaces conceptually a fractured medium with a porous medium, applying equivalent hydraulic parameters. The EPM is commonly used due to its simplicity, as it avoids the detailed description of fractures, but it can only be used if at the representative elementary volume (REV) corresponding with the model grid size (Hassan et al. 2014), there is sufficient fracture density and connectivity (Long et al. 1982).
The HRSs-WLEs, due to their hydrological complexity, require appropriate modelling techniques for reliable simulation. The traditional groundwater models (standalone models) simulate only the saturated zone, applying arbitrary recharge as the driving force input but do not simulate the unsaturated zone water fluxes. Such arbitrary recharge is usually uncertain and highly deterministic when considering the model solution, so also unreliable. The recent new approaches of integrated hydrological models (IHMs) improve that principle. In IHMs, instead of recharge, rainfall minus interception next to potential evapotranspiration are typically used as driving forces, and they are measurable, while the water balances (so also then recharge) are estimated internally by a model based on driving forces and system parameterization.
In the last decades, many IHMs have been developed; each has its own way of adapting the forms of the governing equations for the surface and the subsurface flow systems, the coupling strategy and the solution technique (Maxwell et al. 2014). The most important characteristic of IHMs is that they dynamically couple the surface with the groundwater flow domains across the unsaturated zone, meaning that the model solution of several flow domains is computed within the same single model simulation. The IHMs can be divided, according to the complexity of integration of the surface with groundwater flow domains, into: (1) fully coupled IHMs; and (2) simplified IHMs. The 'fully coupled IHMs' are physically based models in which physically based governing equations of all modelling domains are solved simultaneously. Examples of fully coupled IHMs are CATHY (Camporese et al. 2010), Hydrogeosphere (Brunner and Simmons 2012), MODHSM (Panday and Huyakorn 2004), and PARFLOW (Kollet and Maxwell 2006). All the previously mentioned fully coupled IHMs provide a high level of accuracy and functionality but are also costly in terms of modelling time and intensive in input parameters and computational requirements. The 'simplified IHMs' also dynamically couple the surface with groundwater flow but simplify one or more flow domains' governing equations, being that way more computationally efficient. Examples of simplified IHMs are the MODFLOW-based codes that apply the kinematic wave approximation (KWA; Niswonger et al. 2006) of the onedimensional (1D) Richards equation simulating flow across the unsaturated zone. An example of such 'simplified IHMs' is GSFLOW code, in which the KWA allowed one to avoid highly nonlinear and computationally intensive Richards equation flow simulation. Because of their efficiency, the simplified IHMs are particularly suitable for the regional-scale studies, where the errors associated with the KWA are small, relative to the errors resulting from scaling effects and from a reduced set of parameters representative of a highly complex system (Morway et al. 2012); besides, the simplified IHMs require fewer input parameters. The simplified IHMs can be either: (1) multienvironment or (2) single-environment. A multienvironment IHM involves dynamic coupling of two or more codes such as the coupling of MODFLOW with PRMS in GSFLOW (Markstrom et al. 2008). A single-environment IHM, e.g. MODFLOW 6 ) and MOD-FLOW-OWHM (Hanson et al. 2014;Boyce et al. 2020), can simulate flow across the unsaturated zone by using so-called advanced packages, e.g. UZF (Niswonger et al. 2006), SFR (Niswonger and Prudic 2005), SWR (Hughes et al. 2012), and LAK (Merritt and Konikow 2000), which dynamically integrate the surface system with the groundwater system, across the unsaturated zone.
The latest version of MODFLOW (MODFLOW 6) is an object-oriented framework, which supports the use of multiple models within the same simulation . MODFLOW 6 includes the simultaneous use of MOD-FLOW-NWT (Niswonger et al. 2011), improving simulation of system nonlinearity and MODFLOW-USG (Panday et al. 2013), better handling system heterogeneity within one numerical solution. New operational functions implemented in MODFLOW 6, particularly in the unsaturated zone flow (UZF) and the water mover (MVR) packages, provide new capabilities that handle various deficiencies of the former MODFLOW versions (including those mentioned in the preceding), allowing for more-detailed and more-adequate (compared to the real condition) simulation of complex systems than the previous MODFLOW versions. Finally, these new MODFLOW 6 capabilities allow one to introduce and implement, as in this study, the new cascade-routing and reinfiltration (CRR) concept, which improves streamflow simulation and water balances of MODFLOW 6.
Reinfiltration is a physical process occurring when water from the rejected infiltration (the portion of precipitation that exceeds the infiltration rate) and/or from the groundwater exfiltration (when the water-table rises and exfiltrates to the surface) moves as direct runoff towards nearest stream, lake or any water body. That direct runoff can also reinfiltrate back into the downslope soil or even evaporate before reaching a nearby water body. Although Dogrul and Kadir (2021) state that modelling such a complex process requires detailed information on physical characteristics of soil, land cover, topography, evaporation patterns, etc., the reinfiltration can still be represented in a simplified manner, as proposed in this study.
The research proposed in this study aims to: (1) implement the novel cascade-routing and reinfiltration (CRR) concept in MODFLOW 6, as to the authors' knowledge, it is the only IHM that allows for implementing such concept; (2) design and compute a detailed water balance, applying CRR, in MODFLOW 6; and (3) investigate the dynamics of the surface-water/groundwater interactions in the HRS-WLE, applying new MODFLOW 6 modelling tools.
To address the objectives of this study, the Sardon catchment in Spain was selected because: (1) it is a typical HRS-WLE with challenging surface-water/groundwater interactions; (2) it has the advantage of long-time records of monitored data, which facilitate its on-going research; (3) it has been widely investigated by multiple studies, which allowed for the use of related previous results such as presented by Lubczynski and Gurwin (2005), Mahmoudzadeh et al. (2012), Reyes-Acosta and Lubczynski (2013, Hassan et al. (2014), Francés et al. (2014), Balugani et al. (2017), and Hassan et al. (2017); (4) comparison can be made with the former IHM of Hassan et al. (2014), carried out by the GSFLOW code in the same Sardon catchment; and (5) there is only a little human impact on the land cover and water resources in this area.

Study area
The Sardon catchment (~80 km 2 ) is located in the western part of Spain, ~40 km west of Salamanca city (Fig. 1). The catchment's elevation ranges from 730 m asl at the northern watershed boundary to 860 m asl at the southern boundary. The catchment is composed of weathered and fractured granites with sparse occurrence of other rock types. The southern, western and partly the northern catchment boundaries are marked by outcrops of massive granites with no or very low permeability and granitogneisses with inclusions of shists in the southern part. The eastern boundary extends along a quartzite dyke (Lubczynski and Gurwin 2005). The Sardon River valley is centrally located and oriented approximately in the S-N direction, along the Main Sardon Fault (Fig. 1). The valley has steeper eastern slopes and gentler western ones. The main land use is pasture, as the soil developed from weathered granite is low in organic nutrients, so agriculture activities are rare.
The study area has a semiarid Mediterranean climate, typical for the Central Iberian Peninsula. The precipitation is strongly temporally variable, ranging from ~300 mm year -1 (2009) to >900 mm year -1 (2001). The mean precipitation from 1951 to 2012 was 586 mm year -1 with a standard deviation of 179 mm year -1 (Hassan et al. 2014). The driest months are July and August, with a mean precipitation of <20 mm month -1 , while the wettest months are October and November, with a mean precipitation of >70 mm month -1 . The warmest months are July and August, with a mean temperature of 20 °C and mean potential evapotranspiration of 5 mm day -1 . The coldest months are January and February, with a mean temperature of 5 °C, and the lowest potential evapotranspiration is in December and January, ~0.5 mm day -1 (Lubczynski and Gurwin 2005).
The vegetation of the study area is of savannah type, which in the Mediterranean environment has a form of oak woodland. There are only two types of oak tree species, evergreen oak Quercus ilex (Q.i.) and broad-leafed deciduous oak Quercus pyrenaica (Q.p.), with sparse, ~7% canopy coverage (Reyes-Acosta and Lubczynski 2013). The rest of the area is covered by seasonal grass, occurring only for ~3 months year -1 (typically from April to June), while for the rest of a year, the soil is bare except for some patches of perennial Cytisus scoparius (Scotch broom) shrub.
Different land-cover types affect surface and subsurface hydrological processes, so also affect the catchment system dynamics. That impact must be reflected in the model parameterization; therefore, a land-cover map (with 1-m resolution) was prepared, taking into account two main constraints, surface lithology and vegetation. Following Francés et al. (2014), who mapped surface lithology using two high-resolution multispectral satellite images (QuickBird from August 2009 and WorldView-2 from December 2012), two land-cover classes were defined: rock (mainly granite) outcrops and unconsolidated (mainly weathered) rock deposits (see Fig. 2), referred to as soil. Following Reyes-Acosta and Lubczynski (2013), who classified vegetation using the same two images, three land-cover classes were defined, Q. ilex (Q.i) canopy area, Q. pyrenaica (Q.p.) canopy area, and grass/bare soil. The two maps of Francés et al. (2014) and Reyes-Acosta and Lubczynski (2013) were combined, and as a result, six land-cover classes were obtained as presented in Fig. 2.
The geology is typical for granitic areas, comprising three types of rock materials: weathered granite, fractured granite and massive granite. The thickness of the weathered granite (saprolite), defined by resistivity sounding and boreholes (Lubczynski and Gurwin 2005;Francés et al. 2014), varies from zero at the outcrops of the underlying solid granite The weathered unconsolidated material has a typically sandy fraction with small occurrences of alluvial deposits, which are regarded as the loamy fraction. The deeper the weathering profile, the more gradual transition between weathered and fractured rocks is. The widely outcropping granite shows an abundant network of faults and fractures; this was mapped by Francés et al. (2014), applying a high pass filter on a high-resolution digital terrain model. The Main S-N Sardon Fault ( Fig. 1) controls the course of the Sardon River, while the set of secondary NE-SW faults controls the direction of the tributaries of the Sardon River. The outcrops of massive nonfractured granite are scarce, restricted to local occurrences along the NW catchment boundary.
The hydrogeological framework of the study area was first defined by Lubczynski and Gurwin (2005), where the system schematization consists of three granite layers: (1) spatially discontinuous weathered unconsolidated layer; (2) spatially continuous fissured layer; and (3) massive impermeable bedrock layer. The upper two layers, representing the two aquifers, are hydraulically interconnected. The water table in the study area is unconfined, so its pattern follows the topography, ranging from ~10 m bgs at the elevated areas to less than 3 m bgs in valleys. The shallow water table allowed local farmers to dig artificial ponds supplying water for cattle; some of these ponds dry up in the dry seasons, but those deep enough, with their bottom below the lowest groundwater level, remain water-filled throughout years, clearly indicating the position of the water table. The amplitude of water-table fluctuation is in the order of ~2 m and the water-table fluctuations also influence the intermittent streams controlled by faults and fractures, which drain the aquifers. Towards the dry season, when the water-table declines below the drainage base of a stream, the stream flow ceases; the length of such cessation depends on the amount of rainfall in the preceding wet season and also when the new rains of the subsequent season start.
In the Sardon catchment, there are two automated data acquisition system (ADAS) stations that were implemented to monitor the hydrological variables on an hourly basis (Lubczynski and Gurwin 2005). The first one (Trabadillo) is in the lower, northern part of the catchment, while the other one (Muelledes) is in the upper, southern part, as shown in Fig. 1. The recorded data are climatic variables, particularly the rainfall, air temperature, wind speed, relative humidity, incoming and outgoing solar radiation and barometric pressure. Moreover, there is an automated groundwater monitoring network (Fig. 1), continuously recording groundwater levels hourly since 1994. Additionally, the network includes a flume gaging station, focused on lowflow river measurements (maximum discharge capacity of 145 L s -1 ) that were carried out at the northern catchment outlet in the period of 1997-2001 and afterwards, were extrapolated using water levels recorded in the piezometer at the river channel, near the flume (Hassan et al. 2014).

Conceptual model and water balance components
The first conceptual model- Fig. S1 in electronic supplementary material (ESM)-and the lateral boundaries of the Sardon catchment study area were defined by Lubczynski and Gurwin (2005). The lateral boundaries were delineated along the topographic catchment boundaries ( Fig. 1), which coincide with groundwater divides, except for a small section of the northern boundary, in which the groundwater outflow takes place. Later, Francés et al. (2014) schematized differently the Sardon catchment aquifer geometry, following the general three-dimensional (3D) geological conceptual model of granite aquifers (Dewandel et al. 2006), using a combination of various data sources. They defined two aquifer layers-the saprolite layer and the fissured layer-underlain with the impermeable granite ( Fig. S2 in the ESM). In this study, the conceptual model of Francés et al. (2014) and the lateral boundaries proposed by Lubczynski and Gurwin (2005) were adopted.
The system was conceptualized with two zones, the unsaturated zone and the groundwater (saturated) zone (Fig. 3). Considering the water balance (WB) of the Sardon catchment, the only external water input to the system (source) is the precipitation (P), while the external outputs (sinks) are the evapotranspiration (ET), the Sardon River outflow (q) and the lateral groundwater outflow across the permeable section of the northern boundary (q g ) (Eq. 1).
where ∆S is the change of total catchment storage, equal to the sum of the change of the unsaturated zone storage (∆S u ) and of the change of groundwater zone storage (∆S g ).
The part of P, which falls on the ground surface and is not lost by plant interception (E I ), is called effective precipitation (P e = P -E I ) and in MODLFOW 6 is referred to as 'specified infiltration'. If P e ≤ K v , where K v is the vertical hydraulic conductivity of the unsaturated zone, then the 'net infiltration' (I) is equal to P e . However, if P e > K v , then I = K v and P e -I = RI, where RI is the rejected infiltration (Eq. 2).
The RI, as illustrated in Fig. 3 and presented in Eq. (2), can be (1) evaporated as a sink (RI e ) or (2) transferred downgradient as runoff (RI r ) to be either reinfiltrated back to subsurface (RI i ), or routed to nearby streams or lakes (RI s ). The reinfiltration option has not been available in any MOD-FLOW version yet, so also not in the standard MODFLOW 6 code. This is also the reason why in the current version of MODFLOW 6, (RI i + RI s ) is not estimated separately.
In shallow water-table systems such as the Sardon catchment, when a prolonged and intense rainfall results in a substantial recharge, the water table may rise and exfiltrate to the surface or to a shallow subsurface level, defined by a depth (d surf ) relative to the land surface. The groundwater exfiltration (Exf gw ) in MODFLOW 6 is referred to as 'groundwater seepage' and it is defined as follows: (1) P = ET + q + q g ± △S (2) RI = RI e + RI r = RI e + RI i + RI s (3) Exf gw = Exf e gw + Exf r gw = Exf e gw + Exf i gw + Exf s gw The Exf gw , as illustrated in Fig. 3b, can be (1) evaporated ( Exf e gw ) or (2) transferred down-gradient as runoff ( Exf r gw ) to be either reinfiltrated back to the subsurface ( Exf i gw ) or routed to nearby streams or lakes ( Exf s gw ). As already explained, since the reinfiltration option has not been available in the standard MODFLOW 6 code, the two components Exf i gw + Exf s gw could not be estimated separately in the MODFLOW 6 version used in this study.
The estimate of the ET consists of surface evaporation (E s ) and subsurface evapotranspiration (ET ss ): where E ow is evaporation from open water bodies, in this study assumed as negligible (E ow = 0) where ET u is unsaturated zone evapotranspiration and ET g is groundwater evapotranspiration.
The estimate of q is expressed as follows: where RE s is direct runoff originated from the sum of RI s and Exf s gw , and q B = q gsq sg is baseflow, q gs is groundwater leakage to streams and q sg is streams leakage to groundwater.
The WB of the unsaturated zone is presented as follows: where R g is gross groundwater recharge. The lack of a reinfiltration option in the standard version of MODFLOW 6 and related impossible split of the components RI r = (RI i + RI s ) and Exf r gw = Exf i gw + Exf s gw affected also the WB of the unsaturated zone, as the individual components in brackets in Eq. (8) are inseparable. This was solved by presenting Eq. (8) in the form of three alternative Eqs. (9), (10), (11), depending on the relation between the total reinfiltrated water RE i = RI i + Exf i gw and the rejected infiltration transferred down-gradient RI r = (RI i + RI s ).
(4) E s = E I + RI e + Exf e gw + E ow (5) ET ss = ET u + ET g (6) ET = E s + ET ss = E I + RI e + Exf e gw + ET u + ET g (7) q = RI s + Exf s gw + q B = RE s + q B (8) 10) if RE i < RI r → RI r − RE i = RI r net → P e = ET u + R g + RI e + RI r net ± △S u Note in this study the active infiltration I a = I + RE i net . The WB of the groundwater zone is presented as: The net groundwater recharge (R n ) is defined as:

Driving forces
In this study, the driving forces of the model are effective precipitation and potential evapotranspiration.

Effective precipitation ('specified infiltration')
Precipitation is monitored hourly by two ADAS stations at Trabadillo and Muelledes ( Fig. 1), using ARG100 tipping bucket gauges. As the two records are similar and also because Lubczynski and Gurwin (2005) confirmed similarity between daily records of precipitation at various (more than two) gauge locations within the Sardon catchment, the Trabadillo ADAS station was selected as representative for the daily precipitation of the whole Sardon catchment. The hourly Trabadillo ADAS precipitation records were (11) if RE i = RI r → P e = ET u + R g + RI e ± △S u (12) R g + q sg = Exf gw + ET g + q gs + q g ± △S g (13) R n = R g − Exf gw − ET g corrected for the wind speed following the method outlined by Pollock et al. (2018), aggregated to daily records and assigned as representative for the entire Sardon catchment model.
An experimental tree interception study was earlier carried out in the Sardon catchment by Hassan et al. (2017), who measured interception rates of the two oak tree species, Q.i. and Q.p., in hydrological years 2012 and 2013, and scaled them up to the whole Sardon catchment, applying very high-resolution remote sensing (RS). Furthermore, they extrapolated the tree interception rates in time, using the revised Gash analytical model (Gash et al. 1995;Eq. 15). These tree interception data were used in this study. However, their study did not include the grass interception. Therefore, in this study, the RS-based, revised Gash analytical model was used to derive interception losses of the grass. The revised Gash model assumes that rainfall occurs as a series of discrete events. Each event consists of three periods: (1) wetting up period, when rainfall P is less than the amount of rainfall required to fully saturate the canopy, P' (Eq. 14); (2) saturation period, when rainfall rates are ≥ 0.5 mm h -1 (Gash 1979); and (3) drying out period after rainfall ceases. The rainfall events were defined by a separation of at least 3 h with no rain to allow for complete drying of the canopy.
where E I is canopy interception rate, R is mean rainfall intensity during a day (retrieved from ADASs measurements), S is canopy storage capacity, c is fractional canopy cover, S c = S/c is canopy storage capacity per unit area of canopy cover, ET o is reference evapotranspiration calculated by the FAO Penman-Monteith method (Allen et al. 1998), ET o is mean reference evapotranspiration rate during a day, ET oc = ET o ∕c is mean reference evapotranspiration rate per unit area of canopy cover. In Gash's model, there are two main parameters that are related to the canopy properties, S and c. The leaf area index (LAI) can be used to get c. The LAI~c formula (Steven et al. 1986), used in many hydrological models, e.g. in the HYDRUS-1D model (Šimůnek et al. 2013), was applied (Eq. 16). The LAI is also considered a good predictor of S, as proved by many previous studies (Vegas Galdos et al. 2012;Gómez et al. 2001), in which LAI~S relationships were derived for different crops. In this study, the Menzel (1997) LAI~S formula for a grassland was applied (Eq. 17).
The grass LAI was retrieved from the remote sensing MODIS product (MCD15A3H v061, Myneni et al. 2021), with 4-days temporal resolution and 500-m spatial resolution. The MCD15A3H product includes five layers, one of which is the canopy LAI layer used in this study, further referred to as LAI-MODIS. All the images for the modelled period of 7 hydrological years, i.e. 1 October 2007-30 September 2014, were downloaded, while some gaps due to cloudy images were filled in by interpolation. To select the LAI-MODIS pixel most representative for the grassland of the study area, the LAI-MODIS images (500 × 500 m) were overlain on top of the land-cover map, and the pixel with 99.5% of grass coverage was selected (for the location see dashed-blue rectangle in Fig. 2). To obtain catchment grass-LAI distribution, the temporally variable LAI of the selected grass pixel was retrieved and extrapolated over the "grass/bare soil" class as in Fig. 2.
The daily LAI values were converted to c and S as per Eqs. (16) and (17) respectively, and at each daily time step calculation, next to the daily P (measured by ADASs) and daily ET o (calculated), they were used in Eq. (15) to obtain the daily grass E I .
The final daily interception maps were constructed using the land-cover map (Fig. 2). Out of the six identified classes, for the interception map, four different classes were used, as the classes "Q.i. on soil" and "Q.i. on outcrops" and "Q.p. on soil" and "Q.p. on outcrops" were combined into "Q.i." and "Q.p.", respectively. Furthermore, the land-cover class "outcrops" was assigned zero E I . The final E I for each of the three-remaining land-cover classes (Table 2) were implemented in the numerical model (section 'Hydraulic parameters and boundary conditions').
Considering the assessment of effective precipitation (P ei ), in general, both precipitation (P i ) and interception (E Ii ) may differ per ith land-cover class, while different landcover classes usually have different coverage fractions (a i ) of the total catchment area. Therefore, based on the concept of an area-weighted average, the P ei of a given land-cover class i is calculated as per Eq. (18), while for the whole catchment it is calculated as per Eq. (19): where P ei is effective precipitation of i-th land-cover class; P i is precipitation of i-th land-cover class (in this study, precipitation was assigned as spatially uniform, so P i = P and P e = P − ∑ n i=1 a i × E Ii ; E Ii is the interception of i-th landcover class; and a i is the i-th land-cover-class coverage fraction over the total catchment area for the calculations of P e (note ∑ n i=1 a i = 1 ), i is an index of a land-cover class, and n is the number of land-cover classes.

Potential evapotranspiration (PET)
Potential evapotranspiration is the upper limit of the evapotranspiration from the vegetation canopy and soil that would occur under given energy and infinite water supply. The PET is calculated following Eq. (20) (McMahon et al. 2013;Dassargues 2018), which involves two components: (1) ET oreference evapotranspiration estimated applying the FAO Penman-Monteith model (Allen et al. 1998); and (2) K ccrop coefficient, which can be split into two components, one addressing the evaporation process (K e ) and the other the transpiration process (K cb ). Following the same concept of the area-weighted average as for mapping interception, PET is calculated per each of the six land-cover classes (Eq. 20) and for the whole catchment (Eq. 21).
where PET i is PET of i-th land-cover class, K ei is K e of i-th land-cover class, K cbi is K cb of i-th land-cover class, b i is the i-th land-cover-class coverage fraction over the total catchment area for the calculations of PET (note ∑ n i=1 b i = 1 ), i is an index of a land-cover class, and n is the number of land-cover classes.
The ET o needs meteorological data, particularly the net radiation, wind speed, air temperature and relative humidity (Eq. 6 in Allen et al. 1998). All these data were retrieved hourly from the ADAS station and used to estimate the daily ET o . The ET o was applied as temporally variable but spatially uniform. In contrast, the applied K c was not only temporally but also spatially variable, differing per i land-cover class (Fig. 2).
In the assessment of PET and related definition of a landcover class, if there are trees (or shrubs), it is important to define the area influenced by tree water uptake, that is defined by lateral root extent (LRE). For the "grass/bare soil" (LC1), both K e and K cb were estimated, so K c = K e + K cb , while for the areas influenced by trees within their LREs, whether on soil or on outcrops (LC3-6), only K cb was estimated (K c = K cb ), and for "outcrops" (LC2), K c = K e . The LRE area of Q.p. was assumed to be equal to the ground projection of a canopy, while for the Q.i., the conservative LRE radius multiplier of 2.5 was used following Moreno et al. (2005) who defined Q.i. LRE as 2.5-7 times larger than the canopy radius. Nevertheless, the assumption that the Q.i. LRE is larger than its canopy's ground projection, implied that the areas attributed to these land-cover classes (described by b i coverage fractions) are different from those applied for interception (described by a i coverage fraction).
The K e1 of the "grass/bare soil" was estimated after Balugani et al. (2017), who experimentally defined the evaporation of the "grass/bare soil" (LC1, Fig. 2) land-cover type in the Sardon study area. The daily K e1 of the hydrological year 2010 was obtained by dividing the daily soil evaporation of Balugani et al. (2017) by the corresponding daily ET o and was assumed as valid for all other hydrological years of the model simulation period. The "outcrops" (LC2, Fig. 2) is represented by fractured granites; those fractures are relatively dense and generally open towards the shallow water table, so evaporation takes place but at a lower rate than over the "grass/bare soil" class. Considering the lack of experimental evaporation data on fractured outcrops, the K e2 was arbitrarily approximated as 0.5 K e1 .
To estimate K cb , many studies relate the transpiration rates with vegetation indices such as the normalized difference vegetation index (NDVI) and the soil adjusted vegetation index (SAVI). Gonzalez-Dugo et al. (2009) had derived a linear relationship between K cb and SAVI (Eq. 22), which was later used by Carpintero et al. (2020) for retrieving K cb for Q.i. in the Dehesa region. The same approach was used to get the daily K cb for "grass/bare soil" and for both tree species, Q.i. and Q.p. The SAVI values were retrieved from a series of Landsat-7 TM images (after removing the cloudy images) for the modelled period of 7 hydrological years, i.e. 1 October 2007-30 September 2014, with 16 days temporal resolution and 30-m spatial resolution.
where c = SAVI−SAVI min SAVI max −SAVI min is fractional canopy cover, c max is the c at which K cb is maximal (K cbmax ); K cbmax was assigned as 1 for the "grass/bare soil" and 1.3 for both tree species, Q.i. and Q.p.
To obtain spatial variability of SAVI, the processed Landsat-7 TM SAVI maps were resampled to 1 m resolution and overlain on top of the land-cover map. Due to the very low spatial variability of SAVI, for each land-cover class, an average SAVI value of each of the daily SAVI maps was calculated, and interpolation was performed to get the missing daily values.

Implementation in MODFLOW 6
For setting up the numerical model, the open-source Python script (FloPy) was used (Bakker et al. 2016), which supports different versions of MODFLOW, including MODLFOW 6. The main advantage of FloPy is its flexibility to support all the MODFLOW capabilities, including the packages that are not implemented in the standard graphical user interfaces facilitating MODFLOW model implementation.

Grid design
MODFLOW 6 is based on a generalized control volume finite-difference (CVFD) method, in which a cell can be connected to any number of arbitrary cells ); hence, it supports different unstructured grid types such as triangular nested, rectangular quadtree, rectangular nested and Voronoi grids. The CVFD requirements (Panday et al. 2013;Langevin et al. 2017) are that the line connecting the centers of two cells should: (1) intersect the shared edge at a right angle; and (2) bisect the shared edge. For clarity, the two CVFD requirements are illustrated in Fig. S3 of the ESM. Not all unstructured grid types meet the CVFD connection requirements, but the closer the grid honours the CVFD requirements, the smaller the accuracy losses in the groundwater flow solution. The Voronoi grid (Fig. 4) has the advantage of closely honouring the CVFD requirements (Hesch 2014), particularly when the local grid refinement is smooth (this was achieved by creating a group of transition zones for changing the size of the adjacent cells smoothly, Fig. 4c). Moreover, the ghost node (GNC) package can still be used with the Voronoi grid to minimize the model accuracy losses. Finally, out of all grid types, the Voronoi grid is the most flexible in realistic representation of curvatures of irregular features such as faults or streams (Fig. 4d), so it allows one to optimally represent the spatial system variability.
A combination of FloPy, Python tools and ESRI Arc-GIS was used to build the Voronoi grid (Daoud 2020). For enhancing the accuracy of the model solution, the GNC package was activated; later, it was found that it resulted in ~5% improvement in the model solution. The grid was horizontally discretized in such a way that the smallest cells (width ~15-20 m) were assigned along the streams and became larger further away from the streams (maximum width ~200 m; Fig. 4). The grid had two vertical layers, with thickness defined earlier by Francés et al. (2014) (Fig.  S2 of the ESM). The land surface elevation was retrieved from the 5-m resolution digital elevation model of the Spanish National Center for geographic information (Instituto Geográfico Nacional n.d.). Layer 1 cells, corresponding with outcrops of layer 2, were excluded from the model solution, assigned as vertical through-pass cells (Fig. 4a), using the option "IDOMAIN = -1" in the discretization by vertices (DISV) package.

Hydraulic parameters and boundary conditions
The initial values of the hydraulic and storage parameters are shown in Table 1. The K h = 0.1 m d -1 and anisotropy ratio K v /K h = 1/10 were initially defined in the node property flow (NPF) package as spatially uniform and later adjusted during the model calibration using a group of K-zones for each layer separately. The storage parameters were also preliminarily defined as spatially uniform, i.e. specific storage (S s = 10 -5 m -1 ) and the specific yield (S y = 0.05), in the storage (STO) package and later adjusted during the calibration by using a group of S s and S y zones for each layer separately.
The watershed divide, that surrounds nearly the whole catchment and coincides with the groundwater divide, was assigned as a no-flow boundary, except for the small section (<1 km) at both sides of the Sardon River outlet (Fig. 4), which acts as a lateral groundwater outflow. That outflow  section was assigned as a head-dependent boundary using the drain (DRN) package. The drain elevation, the hydraulic conductivity of the drain's bed (K d ) and the drain's bed thickness (b d ), were assigned as 733 m, 0.05 m day -1 and 0.6 m, respectively. The K d was adjusted during the model calibration as dependent on the K v of the matching cells.
To simulate the flow through the unsaturated zone, the UZF package was activated. The UZF package simulates the vertical unsaturated flow using the KWA of Richards' equation, solved by the method of characteristics (Niswonger et al. 2006). All the active cells of layer 1 and the uppermost cells of layer 2 (those below the outcrop cells) were defined as UZF cells (Fig. 4a).
The UZF parameters are listed in Table 1, where the soil texture-related parameters (θ resid , θ sat , θ ext ) were assigned as spatially uniform based on abundant measurements in the study area (Lubczynski and Gurwin 2005;Hassan et al. 2014;Francés 2015;Balugani et al. 2017). To calculate the actual ET fluxes, the model uses extinction water content (θ ext ) and extinction depth (d ext ) parameters. The θ ext = 0.05 was assigned as spatially uniform and equal to residual water content (θ resid ). The d ext for the land-cover classes "grass/ bare soil" and "outcrops" were assumed as 1 and 0 m respectively, while for both tree species, Q.i. and Q.p. were assigned as 3.7 m based on Canadell et al. (1996). To have one value of d ext per each grid cell, the same area-weighted average concept was used, as presented in Eq. (21). The option for simulating groundwater exfiltration was activated, and the user-specified depth (d surf ) relative to the land surface where the groundwater exfiltration starts, was assumed as d surf = 0.125 m. Note, in MODFLOW 6, the user assigned SURFDEP = 2 × d surf .
The SFR package was activated to simulate the flow interaction between the streams and the groundwater. The streams (the Sardon River and its tributaries) were all defined as SFR reaches in such a way that only one reach was assigned to one grid cell (Fig. 4) with the following data needed: stream's length, width, slope, Manning coefficient, bed level, bed thickness and bed hydraulic conductivity. The streams' length and slopes were calculated using ESRI ArcGIS software, while the stream's width was assumed as 10 m for the Sardon River and 5 m for its tributaries based on the fieldwork observations. The streams' channels were defined as rectangular cross-sections. The Manning coefficient and bed thickness were assumed as 0.035 and 0.2 m, respectively, for all the stream reaches. The bed levels were calculated as equal to the top level of the grid cell that contains the stream reach minus the stream depth. The stream's depth was assumed as 2 m for the Sardon River and 1 m for the Sardon tributaries based on the fieldwork observations. The streams' bed hydraulic conductivity (K b ) values were adjusted during the calibration as being dependent on the K v of the matching cells that contain the reaches.

Implementing the CRR concept
To simulate the direct runoff and to apply the CRR concept, the MODFLOW 6 water mover (MVR) package was used, which has been recently introduced. This package allows moving water from a feature in one package as a provider to a feature in the same package or in another package as a receiver. There are four different options for controlling such water movement . One of these options, "FACTOR", was used in this study to control the transfer of rejected infiltration and groundwater exfiltration (upslope UZF cells' providers), within a given time step, to the adjacent downslope feature(s) (receivers), where that water can be: (1) evaporated ( RI e + Exf e gw ) and/or transferred downslope to be either (2) reinfiltrated back to subsurface ( RE i net ) if the adjacent feature(s) represent UZF cell(s) (not fully saturated yet) and/or (3) discharged as direct runoff (RE s ) into a stream, to contribute to streamflow, if the adjacent feature(s) represent SFR reach(es). However, in MODFLOW 6, there is no automatic way to define the mover fractions (FACTORs) from each UZF provider to adjacent receivers, so the mover fractions were calculated in this study outside MODFLOW 6, applying the multiflow direction (MFD) concept using a combination of Python tools and Microsoft Excel software.
The MFD concept is a raster-based algorithm, which partitions the flow among the downslope neighbouring pixels based on the land surface gradient of the connected pixels (Quinn et al. 1991). The MFD concept is used in many hydrological models and is applied in many GIS applications such as the "flow direction" tool in ESRI ArcGIS software (Qin et al. 2007). Similar concepts as the MFD have already been applied for grid-based models such as the cascade routing tool (CRT) (Henson et al. 2013) of the numerical code GSFLOW (Markstrom et al. 2008). However, in contrast to GSFLOW, the developed CRR concept (Eq. 23; Fig. 4b), allows for surface flow routing of not only regular (as in the CRT) but also of irregular grid cells such as the Voronoi grid applied in this study.
In the CRR concept, the flow fraction (α i,j ) from a cell i to the neighbouring cell j is calculated as follows: where S i,j = (elv i -elv j )/l ij is the slope gradient between the center cell i and j-cells surrounding the cell i (Fig. 4b), elv i and elv j are land surface elevations of cells i and j respectively, l ij is a distance between the center of the cell i and the center of one of j-cells surrounding i-cell, and m is the number of connected j-cells to the cell i; note, the negative value of S i,j means that cell i has a lower elevation than cell j and for such m-connection α i,j = 0 and no flow occurs between that i-j connection. β i,j is a flow partitioning factor used in model calibration, allowing for appropriate partitioning between evaporated water and direct runoff. The α i,j and β i,j range from 0 to 1.
Considering cell i surrounded by m-number of j-cells-so also the m-number of cell connections each described by S i, j -there are three possible flow scenarios: (1) if among all m-connections around an i-cell, only one S i,j is positive (the elevation of the i-cell is higher than the elevation of only one connected j-cell), while all the other S i,j are negative (so they all will be considered as zero in the sum in the denominator of Eq. 23), then α i,j = 1 representing the single flow direction (SFD); (2) if among all m-connections around an i-cell, there are two or more positive S i,j , then the corresponding individual α i,j values will be within the range 0 < α i, j < 1; and (3) if all m-connections around an i-cell are negative S i,j , then that i-cell will be a sink cell where all the water will be evaporated as RI e and Exf e gw , so also contributing to the total surface evaporation (Eq. 4).

Model calibration, validation, and sensitivity analysis
As the first step, a steady-state model was created using averages of 7 years of data, including specified infiltration rates, PET rates, stream stages and observed variables. The calibrated steady-state model was useful for: (1) giving the first indication of the calibration parameters for the transient model; and (2) using its outputs as initial conditions for the spin-up period in the transient model calibration. To initiate the transient model run, which is always a challenge, the following strategy was applied. The daily stresses of the first year of the transient model calibration (1 October 2007-30 September 2008 were duplicated and assigned as an arbitrary spin-up period for the transient model run, consisting of a total period of 8 years (1 arbitrary spin-up year + 7 actual simulated years, from 1 October 2007 to 30 September 2014), equal to 2,922 daily stress periods, each consisting of one time step. The transient calibration was the most time-consuming step of this study-for example, one forward run of the 8-years model took around 14-16 h, using a powerful laptop (Intel Core i7-9 th generation processor and 16-gigabit memory). The transient model was calibrated using the daily values (obtained from the monitoring network described in section 'Study area') of the 14 groundwater head observation points and streamflow at the catchment outlet point (locations are shown in Fig. 1).
The progress of the transient model calibration was controlled by applying graphical and statistical comparisons between the simulated and the observed values to minimize the differences in 14 head points to achieve the RMSE ≤ 1 m set as a calibration target. The calibration also involved the comparison of the simulated with the flume-estimated catchment outflows using the volumetric efficiency (VE), a metric for flow observations ranging from 0 for a poor fit to 1 for a good fit. The VE was formulated from the Nash-Sutcliffe coefficient (NSE) to overcome the NSE failure to represent useful evaluation when NSE < 0 (Criss and Winston 2008). The model was considered sufficiently calibrated if VE ≥ 0.5.
A sensitivity analysis was performed at the early stage of the model calibration on most of the model parameters.
Hereafter, the sensitivity of the two parameters affecting the CRR concept is presented: (1) K v , as it controls the rejection of the specified infiltration; and (2) the flow partitioning factor (β i,j ), as it directly controls the RI + Exf gw partitioning into ET and q. Table 2 presents the catchment's yearly interception contributions of different land-cover classes and total catchment yearly interception, presented as % of yearly precipitation (P). The largest contribution of 'grass/bare soil' land-cover class, which has a low interception rate, is because of its large spatial coverage of the catchment area.

Driving forces
Considering PET, Fig. 5 presents average monthly variability of K c per land-cover class for one hydrological year, i.e. from 1 October 2009 to 30 September 2010, while the average monthly K e and K cb , used in the calculation of K c , are presented in Table S2 of the ESM. The largest K c is observed at all land-cover classes in April, with the highest K c = 1.28 for the evergreen Q.i. LC3 (Fig. 2). Expectedly, in the first part of the hydrological year, the K c of evergreen Q.i. was larger than of deciduous Q.p., but in the dry season, the K c of the deciduous Q.p. exceeded the Q.i. The lowest K c was for the outcrop LC2 class. Table 1 shows the calibrated transient parameters next to the initially assigned parameters. A number of parameters have not been changed during the transient model calibration, while others such as K h , K v , S y and β i,j , were adjusted; the zonal distribution of K h , S y , and S s for both layers is shown in Fig. 6, whereas β i,j remained spatially uniform.

Calibration results
The graphical comparison between the simulated and observed heads at the 14 observation points during the entire simulation period is presented in Fig. S4 of the ESM. Some piezometers such as PPNO and PGTMO show a good match between the observed and the simulated groundwater heads, while others such as PMU1, show a relatively poor match due to a time shift between the simulated and the observed heads. Some other piezometers such as PGBO and PGJO, show faster responses of the simulated than the observed heads; however, in general, there is an acceptable match (RMSE ≤ 1 m) between the observed heads and the simulated heads. The overall RMSE of the entire model is 0.67 m, while the individual RMSE is in the range from 0.2 to 0.98 m.
The data of the Sardon River flow at the catchment outlet (Fig. S5 of the ESM) was useful in controlling the main river outflow in dry seasons, although uncertain because: (1) it was restricted to low-flow conditions of ≤ 0.145 m 3 s -1 due to limited flume capacity; therefore, a wet season match between the flume discharge estimate (black line) and the total river outflow (red line) was impossible, except for short recession periods with discharges declining below 0.145 m 3 s -1 , as for example in April 2009; (2) in the simulated years 2008-2014, the available stream discharge data were not directly measured, but derived from the relation (Hassan et al. 2014) defined prior to 2008, between the flume discharge and the automated water level measurement in the piezometer located in the river channel near the flume, indicating either surface water or groundwater levels, depending on whether surface channel flow was present or not; and (3) there is an unknown quality associated with the flume maintenance carried out by the local farmer. An additional complication in calibrating the model using stream discharges Table 2 Yearly interception (E Ii ) presented in % of precipitation (P) per land-cover class (Grass/bare soil, Q.i. and Q.p.) and for the whole catchment; a i is the land-cover-class coverage fraction over the whole catchment area that is applied for the interception estimate, where i is index of land-cover class (LC) as presented in Fig. 2;E a i I i = E Ii × a i is a contribution of a given land-cover-class interception to the total catchment interception; note: E I3 = E I4 = E I3, 4 , and E I5 = E I6 = E I5, 6 ; also a 3 + a 4 = a 3, 4 and a 5 + a 6 = a 5, 6 ; besides, "outcrops" E I2 = 0, so despite a 2 = 0.22,E a 2 I 2 = 0 , and therefore it is not listed in Table 2 is that the Sardon River follows the Main Sardon fault, and when the surface-water flow ceases in the river channel, the outflow in the subsurface continues along the Main Sardon fault zone underlying the surface river channel, as also confirmed by tracer tests (Lubczynski and Gurwin 2005). That simulated outflow is reflected in the model by baseflow (green line). In that respect, it is remarkable that in every simulated year, wet or dry, that baseflow rate was similar, with minima in peak dry seasons not lower than ~0.03 m 3 s -1 (= ~30 L s -1 ). Therefore also, the exact match between the flume discharge estimate (black line) and the simulated-by-SFR-package baseflow (green line) was not expected; despite that, as the subsurface channel outflow was relatively low, the VE between the simulated streamflows and flume estimates was 0.48, so reasonable.

Water balance
The daily water balance (WB) of the model simulation was exported from the MODFLOW 6 output files and aggregated into yearly rates, as presented in Table 3, and also averaged over the entire period of the model simulation as illustrated in Fig. 7. Considering 7 year WB means, the only input (Eq. 1) is P = 523.4 mm year -1 , while the outputs are: ET = 72.0% of P, q = 24.7% of P, and negligible q g = 0.2% of P. The ET was substantially larger than q in every hydrological year (Table 3), with an average q/ET ratio of 0.34 over the entire simulation period, largely due to the intermittent character of the Sardon stream network. The ET consists of five components (Eq. 6), three surface components (E s ), which totally contribute 28.7% of ET (E I = 12.5%, RI e = 14.6%, and Exf e gw = 1.6%), and two subsurface components (ET ss ), which contribute 71.3% of ET (ET u = 60.1%, and ET g = 11.2%). The q consists of three components (Eq. 7), the q B = 13.6% of q and RI s + Exf s gw together representing RE s = 86.4% of q, as their separate estimate is not possible in the current version of MODFLOW 6.
In the unsaturated zone (Eq. 8; Fig. 7), the main input is P e (91.0% of P) and the two (nonseparable in the current version of MODFLOW 6) water input components ( RI i + Exf i gw ) = RE i net jointly represent 8.4% of P. The unsat-(a) ( b) 0 1 2 3 0.5 km Fig. 6 Spatial distribution of calibrated a K h and b S y and S s for both layers urated zone outputs R g , ET u and RI were 29.2, 43.2 and 22.3% of P respectively. The net infiltration ( I = P e − RI ) to the unsaturated zone, was 68.7% of P (Fig. 7). However, when calculated by applying the CRR concept, I increased by RE i net (8.4% of P), so the final active infiltration ( I a = I + RE i net ) was 77.1% of P. In the groundwater zone (Eq. 12; Fig. 7), the main input is R g (29.2% of P) in addition to the negligible water input from the streams q sg (0.8% of R g ). The main output from the groundwater zone is Exf gw (65.3% of R g ), which emphasizes the significance of groundwater exfiltration in the Sardon catchment. The groundwater loss to the streams, q gs (12.3% of R g ) was stable along the simulated years (Table 3).
The R n showed specific aquifer dynamics (recharge/discharge conditions) in response to different meteorological conditions (dry/wet years). The dry years with relatively low P, such as 2009 and 2012, showed negative yearly mean R n , while the wet years with high P, such as 2010 and 2014, showed positive yearly mean R n (Table 3), although even in such wet years, that R n was only 11% of R g .

Spatial distribution of water fluxes
The spatial distribution of the selected groundwater fluxes (R g , Exf gw , ET g , and R n ) was extracted from the model output for two contrasting hydrological years: 2009 (dry year) and 2010 (wet year), as presented in Fig. 8. For both years, the R g was high (>100 mm year -1 ), not only in the flat elevated recharge areas, but also in the drainage areas adjacent to the Sardon River and its tributaries, typically aligned with faults adjacent to high fracture density zones, where the water table was shallow. In all these areas, there are not only groundwater recharge (R g ) but also groundwater discharge processes (ET g , Exf gw ), so only the R n (Eq. 13) shows the true, net water input into the aquifer.
The R g is more dependent on local conditions (e.g. soil permeability, interception, etc.) than on the position in the catchment flow system. The Exf gw and the ET g follow the same pattern as the R g , because both are primarily dependent, in MODFLOW, on water-table depth, with the highest values in the drainage areas and in the flat plateaus, where the water table is the shallowest. The spatial distribution of R n is dependent on the summative action of R g , Exf gw and ET g (Eq. 13). The R n was generally positive (recharge areas) in the NW part of the catchment, even in the dry year 2009, but negative (discharge areas) in the southern part, even in the wet year 2010. It is remarkable that the R n in most of the drainage network lines was either positive or negative, depending not only on the position of a section of a stream in the hydrological system but also on the year analyzed; obviously, in the wet year 2010, the R n was more positive than in the dry year 2009. Figure 9a,b shows the temporal variability of groundwater fluxes and evapotranspiration fluxes, respectively, both during the dry year 2009 and the wet year 2010. It can be seen that Exf gw is correlated with R g . The maximum daily values of R g and Exf gw were in winters (December-February) of each year when P was high. In contrast, during summers (June-September), R g and Exf gw were very low due to very low or scattered P (Fig. 9a), but ET g was active and had its maxima in springs (March-May) when the water table was are grouped as RE i net (total reinfiltrated water) while the two surrounded by blue dash-line rectangles also do not have separate values and are grouped as RE s (direct runoff); b box plot showing the ranges of each WB component; black circles indicate outliers out of the extent of the whiskers the highest, while minima in autumns (October-November) when the water table was the lowest. The R n (Eq. 13) had positive peaks (recharge conditions) during winters (December-February), when R g and Exf gw were high, but negative (discharge conditions) in springs (March-May), when R g and Exf gw were moderate but ET g was relatively high. During summers (June-August), the R n was zero or very small negative, mainly due to negligible R g , so also negligible Exf gw and small but relevant ET g .

Temporal variability of water fluxes
It is remarkable that the temporal distribution of ET u was different from that of ET g . ET u was not dependant on the water table as it was the case for ET g but was mainly dependent on PET and rainfall occurrence. The highest ET u (Fig. 9b) was in late springs (April-June) when the soil moisture was high, with peaks of ~1.8 mm day -1 in dry years and ~2.5 mm day -1 in wet years and the lowest (0.05 mm day -1 ) in late summers (August-October) when rain was scarce, so low soil moisture, and in winters (December-February), when PET was low. The temporal variability of RI e and E I followed the temporal variability of P in each year (Fig. 9b), which was expected as both fluxes originated from P, so the correlation with P was forced. The largest RI e and E I were observed in the late autumn and winter periods (November-February), while the smallest in summer (June-September) because of negligible P. The temporal variability of the total ET followed the daily variability of its largest component, i.e. ET u , except for winter (December-February) months, when RI e was dominant.
The temporal streamflow variability of the Sardon River at the catchment outlet is presented in Fig. S5 of the ESM. The characteristic feature of the outflow pattern is that there is a fast river response to rainfall due to the shallow water table, the relatively high permeability of the soil and dense fracture system, matching the dense drainage network of streams. These factors, in addition to low catchment storage, imply a short groundwater travel path and a fast drainage of the aquifer system. Another interesting observation is that the direct runoff (RE s ) starts every year around October, when heavy rainfalls start, and remains very active till February when it begins to decline, to become negligible in the summer months. In the winter months, the rains fill up the storage of the catchment, which then often acts as an overflowing reservoir. Finally, while comparing simulated baseflow (q B ) of dry year 2009 with wet year 2010, it can be noticed (Fig. S5 of the ESM) that the dry season baseflows in both years were the same (~0.03 m 3 s -1 or 30 L s -1 ) while the wet year peak baseflows slightly differed from ~0.05 m 3 s -1 in dry 2009 to 0.07 m 3 sec -1 in wet 2010. Considering that in dry summers, the surface river flow was not present, that 30 L s -1 baseflow is attributed to subsurface catchment outflow through the Main Sardon fault zone underlying the Sardon River at the outlet point, as also notified by Lubczynski and Gurwin (2005).

Sensitivity analysis
The 10-time increase in K v decreased RI r net by 88% and increased Exf r gw by 34%, meant that in total RI r net + Exf r gw increased by 12%, leading to an increase in RE i net by 57% (i.e. from 8.4% of P to 13.2% of P). In contrast, the 10-time decrease in K v increased RI r net by 340% and decreased Exf r gw by 95%, meant that in total RI r net + Exf r gw decreased by 19%, leading to a decrease in RE i net by 97% (i.e. from 8.4 to 0.5% of P).
The changes in the flow partitioning factor (β i,j ) showed a significant effect on the partitioning of the two surface WB Fig. 9 Daily variability of different water fluxes: a groundwater zone fluxes; b evapotranspiration fluxes over the two contrasting hydrological years, dry 2009 and wet 2010; a hydrological year starts from 1 October of the previous year and ends on 30 September of the specified year components, ET and q. The increase in β i,j by 0.2, i.e. from 0.8 to 1, where all the rejected and exfiltrated water is routed to streams, decreased ET by 11% and increased q by 42%. That high q (q = 31.5% of P) was not realistic for such a water-limited environment and also showed poor streamflow calibration (VE < 0.5). In contrast, the decrease in β i,j by 0.2, i.e. from 0.8 to 0.6, increased ET by 11% and decreased q by 43%. That low q (q = 11.3% of P) showed acceptable streamflow calibration but the ET in some years exceeded the PET, which eliminated that option.

Discussion
The Sardon catchment is an area particularly suitable to investigate dynamic processes of surface-water/groundwater interactions so also to test related novel concepts and tools such as those applied in this study because of: (1) an excellent database with long time-series automated measurements; (2) hydrological complexity with distinct and fast surface-water/groundwater interactions due to: (a) shallow water table, low storage and high hard-rock permeability; (b) Mediterranean water-limited environment (WLE) with clear wet (with intense showers) and dry (with droughts) seasons and savannah oak woodland interacting with groundwater; and (c) dense drainage network, hydraulically connected with the aquifer system; and (3) remoteness of the study area, implying negligible human impact on hydrological processes.
The dynamics of the surface-water/groundwater interaction in the Sardon catchment was earlier simulated through numerical models by Lubczynski and Gurwin (2005) and Hassan et al. (2014). The study described here also attempts to simulate surface-water/groundwater interactions but it: (1) improves the driving forces input; (2) uses a more realistic conceptual model (after Francés et al. 2014); and more importantly, (3) implements, in MODFLOW 6, a novel CRR concept of surface-water reinfiltration, allowing for more detailed and more realistic water balance estimates.
Driving forces of IHMs, together with model parameters, decide about model performance. P, E I and PET are the water fluxes typically used in IHMs, so they were also used in this modelling study, as driving forces. The P was assigned as spatially uniform, based on automated rainfall measurements, and improved by undercatch correction according to Pollock et al. (2018). That uniformity was justified by the relatively small size of the study area and small elevation differences, but also by earlier investigations (Lubczynski and Gurwin 2005), which confirmed low P spatial variability in that catchment.
Plant interception (E I ) is often underestimated or even neglected by modellers, but it can represent a significant percentage of P-for example, in densely forested areas, interception can exceed 20% of P (Van Stan et al. 2016;Grunicke et al. 2020). In the Sardon study area, characterized by savannah land-cover type, with only ~7% oak tree canopy cover (Reyes-Acosta and Lubczynski 2013 and temporary grassland, active-green only for ~3 months year -1 , the E I was low (from 8.9 to 12.4% of P, Table 2), although still relevant in WB. The mean E I (2008-2014) was ~9%, so relatively low, as for the Dehesa (in Spanish) or Montado (in Portuguese) land cover (Pereira et al. 2009) to which the Sardon study area belongs. The relatively low E I in the Sardon catchment was due to (Table 2): (1) low tree density (a i = 0.07), with larger contribution of low E I of deciduous Q.p. (a i = 0.05) than high E I of evergreen Q.i. (a i = 0.02); (2) dominant catchment contribution of 'grass/bare soil' land-cover type (a i = 0.71) with low E I ; (3) and large catchment contribution of 'outcrops' land-cover type (a i = 0.22) with E I = 0.
The E I of oak trees was estimated in this study following the experimental study of Hassan et al. (2017) in the same Sardon catchment, which however, did not include the grass E I . In this study, the grass E I was estimated based on Gash's model (Gash et al. 1995), where c and S were defined from their relations with LAI, implemented by remote sensing (RS). That grass E I estimate involves uncertainty, as the Menzel S~LAI relation was defined in a bit different climatic and soil conditions than the Sardon catchment. However, even with such a simplifying assumption, the proposed approach is valuable, as it involves realistic temporal E I variability constrained by true LAI temporal variability (presented in Fig. S6 of the ESM), which anyway is better than an arbitrary grass E I estimate.
Another driving force of an IHM is potential evapotranspiration (PET), which represents the maximum ability of a system to evapotranspire water under unlimited availability of water with the given energy supply. Unfortunately, in the scientific world, there is still no consensus about the way how PET should be estimated. In this study, PET was defined from ET o × K c (McMahon et al. 2013;Dassargues 2018), assuming K c = K e + K cb (Allen et al. 1998), where respective evaporation and transpiration K c -components were land-cover dependent (Fig. 2) and temporally variable (Fig. 5). For ET o there is a straightforward calculation protocol (Allen et al. 1998), while for K e and K cb there are no specific guidelines for the natural environment; guidelines are only available for agricultural crops. As such, the K e was assigned following experimental evaporation measurements in the Sardon catchment by Balugani et al. (2017), while the K cb was approximated by the linear relation with the SAVI index, following Gonzalez-Dugo et al. (2009) and both were spatially extrapolated.
The spatial PET variability, similarly to E I , is dependent on the land-cover type, in particular on density of trees, which are the largest water consumers, and on species type. In the Sardon study area, the spatial coverage of oak trees was very low, i.e. 13 stems per ha with only 7% of canopy cover (Hassan et al. 2017), so oak trees present a low contribution to PET. However, the oak trees, despite limited density, may locally play an important discharging role. Both Q.i. and Q.p. are phreatophytic trees, adapted to harsh drought conditions through hydraulic redistribution skills (Lubczynski 2009;Reyes-Acosta and Lubczynski 2014). Such trees have typically two types of root systems: (1) shallow lateral roots, focusing on shallow soil moisture uptake; and (2) deep taproots, focusing on groundwater or its capillary fringe. If shallow soil moisture is available, phreatophytes use it, but once soil moisture is depleted, e.g. during droughts, phreatophytes switch from shallow lateral root water uptake to groundwater uptake by deep taproots. Many tree species, particularly in arid and semiarid areas, are often exposed to long droughts, developing adaptation skills to survive, which however, create various challenges in the quantitative assessment of the hydraulic redistribution (Lubczynski 2009). One of the largest challenges in the spatio-temporal assessment of PET, so also in the water uptake by trees, in the WLE, is attributed to the size of an area influenced by tree water uptake; typically approximated by lateral root extent (LRE). This is because: (1) within LRE area extents, trees typically capture and transpire more water than in adjacent areas (e.g. grasslands); (2) different tree species take up water from the subsurface with different temporally variable rates; and (3) different tree species have different LREs, ranging from smaller than the ground projection of canopy area, to much larger.
The experimental LRE assessment is cumbersome. The LRE can be approximated by geophysical imaging methods, e.g. by electrical resistivity and/or radar imaging (Attia al Hagrey 2007), but its accurate definition is still only available through root excavations. In this study, none of such assessments was possible to carry out; therefore, the LREs of Q.i. and Q.p. oaks were assigned based on literature root information about these two tree species in similar Dehesa/Montado land-cover areas, applying a fixed proportion between the radiuses of LRE and the tree canopy projection (Schenk and Jackson 2002;Hardiman et al. 2017). The LRE was defined automatically for all tree species of the Sardon catchment applying very high-resolution digital image processing. The application of literature-based proportion values between the LRE and tree canopy radiuses, even for the same species but in different study area environments of growth, brings uncertainty that can be mitigated only by dedicated, experimental root studies.
The model calibration was performed using the manual trial and error method, which in contrast to automated parametric optimization (e.g. with PEST or UCODE), allows for reasonably fast processing of numerically complex models such as the Sardon IHM on a standard PC. Also, the zonally uniform distribution of hydraulic parameters (K, S), particularly appropriate for block-faulted, granitic hard-rock areas (Francés et al. 2014;Francés 2015), helped to keep the simulation to a reasonable ~16 h run time, i.e. 3-4 times faster than applying automated parameter optimization code. That zonal parameter distribution also allowed one to have control over the parameter adjustment, often intuitive based on general hydrogeological knowledge about the study area modelled.
The steady-state model calibration provided approximate estimates of unsaturated and saturated zone parameters, but the calibrated state variables were unusable as initial conditions of the transient model. The same problem was reported by Hassan et al. (2014), El-Zehairy et al. (2018, and Lekula and Lubczynski (2019). To initialize the transient model calibration, a simple and efficient spin-up solution was proposed, in which the first year of available data was used to warm up the model, but afterwards, the same year was used in the true-transient simulation. In that way, no data were 'sacrificed' as a spin-up period.
The transient calibration showed a good match between the observed and the simulated variables with RMSE ≤ 1 m for groundwater heads and VE ~ 0.5 for streamflow. Considering the calibration discrepancies of the groundwater heads, they could be due to: (1) errors in model conceptualization; (2) errors in model parameterization; and (3) the effect of grid-scale head variability related to grid-scale heterogeneity.
Considering the calibration discrepancies of the streamflow at the catchment outlet, the good match was not expected because: (1) the flume measured only outflows up to 145 L s -1 , so there was a minimal spectrum of discharges potentially to be used for calibration purposes; (2) the flume measured only the surface outflow, while part of that outflow corresponds to the subsurface outflow along S-N Main Sardon fault zone matching the Sardon River (Lubczynski and Gurwin 2005); that outflow was not recorded by the flume in the dry season when the surface river flow ceased, but the model simulation indicated a stable record of outflow of ~30 L s -1 , regardless of whether in the dry or wet year; (3) there was inaccuracy in the flume discharge measurements and discharge extrapolation; and (4) eventual errors mentioned in the preceding paragraph also as 1, 2, and 3.
The CRR concept was able to simulate the reinfiltration process in a simplified manner, based on: (1) the flow directions, determined by the area's elevations and slopes (Eq. 23); (2) the saturated vertical hydraulic conductivity (K sat ) of the unsaturated zone (herein K sat = K v ); and (3) a flow partitioning factor (β i,j )RI e + Exf e gw . The CRR includes three water components; one of them, the reinfiltrated water, is controlled by the K v , while the partitioning of the other two, evaporated water ( RI e + Exf e gw ) and direct runoff (RE s ), is controlled by the β i,j . The sensitivity analysis of β i,j showed its significant impact (section 'Sensitivity analysis') on ET and q. In this study, the β i, j control was limited to the low flows at the catchment outlet. In follow up studies, it is recommended to have more constrain on β i,j , by for example total streamflow, evapotranspiration, or soil-moisture estimates. The implementation of the CRR concept in MODFLOW 6 substantially affected the WB of all model domains. In the standard MODFLOW 6 (i.e. without CRR), the rejected infiltration (RI) and/or groundwater exfiltration (Exf gw ) were treated as sinks (i.e. as evaporated water), but in reality, that water can still reinfiltrate back to the unsaturated zone (Eq. 8) or move as runoff to be discharged in nearby water bodies. The CRR implementation in the Sardon MOD-FLOW 6 model changed its WB as follows: (1) ET from 93.2 to 72% of P; (2) q from 2.9 to 24.7% of P; (3) R g from 19.9 to 29.2% of P, (4) Exf gw from 11.4 to 19.1% of P. It is remarkable that the ET simulated with CRR (in this study) matched well the ET measurements carried out with the eddy covariance tower by Balugani et al. (2017) in the same Sardon catchment study area; moreover, the WB components obtained with MODFLOW 6 with CRR were much closer to the GSFLOW- Hassan et al. (2014) WB components (73% of P and q of 27%) than without CRR. In both model solutions, the ET was the largest discharging WB component, but was also relatively low as for the WLEs, where P is typically comparable with ET, and q is small or negligible. For example, in the HRS-WLE Hout catchment in South Africa, simulated by MODFLOW-OWHM (Ebrahim et al. 2019), the ET was ~92% of P, likely because of the drier climate, active agriculture activity, deeper water table and much sparser stream network (outflow only 3.3% of P) than in this study. In the extreme case of still drier conditions and a deeper water table (>60 m) of the Central Kalahari Basin study (Lekula and Lubczynski 2019), the yearly ET was approximately equal to the yearly P.
One of the most important characteristics of the Sardon catchment is the large contribution of Exf gw in the WB. The importance of Exf gw is particularly distinct in areas with a shallow water table and stormy rainfall, implying fast water-table rises towards the surface where groundwater can exfiltrate. The Exf gw events are closely related to R g events, while the R g -Exf gw difference, also referred to as effective recharge (Markstrom et al. 2008), together with ET g , constrain the R n (Eq. 13), which represents the critically important WB component, determining either replenishment or decline of groundwater resources. The implementation of CRR increased the contribution of R g in the WB, but also increased the fraction of Exf gw in R g from 57.4 to 65.3%, so it became almost the same as 69% in the GSFLOW study by Hassan et al. (2014). In other available IHM studies (only IHMs are able to simulate Exf gw ), that fraction was typically lower. For example, in a GSFLOW modelling study in a humid (P = ~1,250 mm year -1 ) Miho hard-rock catchment in South Korea (Joo et al. 2018), the fraction of Exf gw in R g was 'only' 59%, mainly because of the deeper water table, while in the extreme case of the Central Kalahari Basin with a very deep water table, the Exf gw was negligible (Lekula and Lubczynski 2019), so the R n could be estimated from R g -ET g . Considering the relevance of Exf gw in WB, particularly its impact on the R n , and also very little scientific information about it, more research, particularly experimental, needs to be dedicated to Exf gw to better understand its behavior in various hydro(geo)logical conditions and to quantify it.
The spatial distribution of R n showed a mosaic pattern, where the recharge and discharge areas were close to each other (Fig. 8). Such a pattern is characteristic for HRSs, as also acknowledged by Hassan et al. (2014). The large spatial variability of the groundwater fluxes in the Sardon catchment (in both dry and wet years, Fig. 8) is mainly because of: (1) hard-rock, granite tectonics of the study area, with faults and fractures, but also spatially interchanging granite outcrops with saprolite of variable thickness; (2) hilly topography enhancing the short travel path and related short residence time of groundwater; (3) shallow water table; (4) dense drainage network; and (5) spatially variable savannah land-cover type with sparse phreatophytic trees locally uptaking groundwater.
The semiarid climate of the study area, with distinct dry and wet seasons, low catchment storage, shallow water table and fast replenishment due to relatively high permeability, imply large temporal variability of nearly all water fluxes. Obviously, that variability is constrained by large daily and seasonal P variability, i.e. scarce or no P in dry summer seasons and highly variable P in the wet winter seasons, with P peaks enhancing maxima of R g , Exf gw , RI and q. The temporal patterns of ET and ET ss are similar because in the Sardon study area, the ET ss is the main contributor of ET. It is remarkable and characteristic for the Sardon catchment that the ET and ET ss peaks occur in early spring, i.e. in April-May (Fig. 9b), before the solar radiation maxima. Such early peaks of ET are due to substantial soil moisture availability remaining still after the wet season and rapidly increasing PET in spring towards summer. Another interesting ET-related observation, is that there are some scattered ET peaks also in the winter months, which are related to flooding episodes when RI e is the main ET contributor.
The temporal variability of ET and ET ss in the Sardon catchment is similar to other savannah study areas. For example, Campos et al. (2013) investigated the evapotranspiration (ET EC ) of the evergreen Q. ilex oak woodland (Caceres, southwestern Spain) in the period of 2004-2008, applying an eddy covariance tower. The ET EC peaks were in May and were similar to the ET rates reported in this study. Another example is the study of Miller et al. (2010) in a semiarid California oak savanna site, characterized by two different vegetation types, grass/soil and deciduous blue oak trees (Q. douglasii; Q.d.), where the highest ET occurred also in spring (March-May) when both the trees and grasses were active. That study emphasized the importance of vegetation not only as a contributor to ET but also to ET g . Miller et al. (2010) showed that in the wet year (2006) the grass/soil ET was 61.0% of the total ET and the Q.d. ET g was 39.0% of the total ET, while in the dry year (2008), the grass/soil ET was 32.0% of the total ET but the Q.d. ET g was 68.0% of the total ET.
The importance of WLE vegetation in ET g formation has been extensively addressed by Lubczynski (2000Lubczynski ( , 2009, who pointed out that the ET g is primarily dependent on the energy supply. However, in this study, the simulated ET g maxima (Fig. 9b) were already in spring when the water table was the highest, but negligible when the water table declined to the lowest position at the end of the dry season. Such an ET g pattern associated with MODFLOW 6 (but also with any previous MODFLOW versions) primarily depends on water-table position, and is likely unrealistic because: (1) in dry seasons, when soil moisture is depleted and grass withered, the 'dry soil layer' is developed in the upper part of a soil, which allows for substantial evaporation by vapor flow; that rate is much larger than predicted by models commonly used in hydrology (so also by MODFLOW 6) and strongly depends on PET and less on water-table depth (Balugani et al. 2021), therefore the evaporation being the highest in the peak dry season; (2) phreatophyte trees, such as Q.i. and Q.p., have access to groundwater through deep taproots, so they are primarily dependent on solar energy supply, which is also the largest in the peak dry season (Reyes-Acosta and Lubczynski 2014). Both observations suggest that the largest ET g in the Sardon study area is in the peak dry season (July/August), when however, the ET g simulated by MODFLOW 6 is relatively low, just because of its primary dependence on the declining water table. More appropriate behavior of ET g was reproduced by Francés (2015), who applied the MARMITES-MODFLOW (MM-MF) model in the La Mata sub-catchment of the Sardon catchment. In that model, the ET g was primarily dependent on PET, by first using the available soil moisture of the unsaturated zone and then completing PET requirement through the ET g under the condition that roots are deeper than the water-table depth. The MM-MF solution exhibited the highest ET g in the dry season between July and October, when the ET u was low, which was in agreement with the observations of Reyes-Acosta and  and Reyes-Acosta (2015). Considering the primary dependence of MODFLOW 6 upon the water-table depth, and related underestimation of the dry seasons' ET g , particularly distinct in dry WLEs, further studies on the conceptualization and implementation of ET g in MODFLOW 6 are recommended.
MODFLOW 6, thanks to CVFD, has the advantage of using any kind of structured or unstructured grids. In this study, the most flexible, Voronoi (Fig. 4) grid was applied. Such grid has optimal ability to realistically represent the important hydrogeological features, such as the curvatures of streams, lakes, seashore, faults and sharp boundaries of block heterogeneities, observation points etc. For example, the Sardon streams were represented by the minimum grid cells' width (~15-20 m), nearly close to the real width of those streams (~10 m), which enhanced the quality and accuracy of the simulation. Another advantage of the Voronoi grid is that it honours the CVFD connection requirements, which reduces the need for corrections (errors in simulated heads and flow due to violation of the CVFD) using the GNC package. That reduction was confirmed in this study by activating and deactivating the GNC package, where the GNC deactivation showed a low effect on the model solution (3% change in the overall groundwater heads' RMSE).
MODFLOW 6 introduces new concepts that also improves UZF package performance; the three main improvements are as follows: (1) the UZF package assignment is no longer restricted to the uppermost layer; (2) the former approximation of the residual water content (θ resid ) by specific retention (S r ), which in reality are two different parameters, is now fixed so that the two can be represented by their adequate, true estimates; and (3) the groundwater exfiltration (Exf gw ) has a more flexible representation because of the newly introduced d surf option, allowing the Exf gw to start from below the land surface.
The application of CRR in MODFLOW 6 introduced new WB components, RI e and Exf e gw , as part of the total E s (Eq. 4). That E s , and so also E I (Eq. 4), is ignored in the simulation of ET, as the UZF package of MODFLOW 6 simulates only ET ss (Eq. 5), assuming that ET = ET ss . The problem of that assumption is that the total ET, after adding the E s to the simulated ET ss (Eq. 6), may exceed PET. It is therefore recommended to consider E s as part of the total simulated ET in the coming versions of MODFLOW 6.
The SFR package simulates 1D surface flow with no stream storage, which could lead to simulation errors if the water travel time of a stream network (from the starting point of the network to the outlet) is much greater than the model time step. Besides, the current MODFLOW 6 version is restricted to only rectangular channel cross-sections. Those limitations are handled in the Surface-Water Routing Process (SWR1; Hughes et al. 2012), incorporated in MOD-FLOW-OWHM (Hanson et al. 2014;Boyce et al. 2020), which simulates 1D and two-dimensional (2D) surface-water flow and has various options of channel cross-sections. The SWR1 capabilities (or similar) would improve the streamflow simulation if linked to any of the next versions of MOD-FLOW 6.
The very useful MVR package is introduced in MOD-FLOW 6. That package is designed to move water from a feature in one package to a feature in the same or another package. In this study, it allowed implementing the new CRR concept, in which RI and Exf gw can be: (1) evaporated and/or moved down-gradient to the adjacent feature(s) to be either (2) reinfiltrated back to the subsurface if adjacent feature(s) represent UZF cell(s); and/or (3) discharged in a nearby stream, if the adjacent feature(s) is an SFR reach(es) representing a stream or other sink-water-body (e.g. lake or sea). However, the disadvantage of the current version of MODFLOW 6 is that splitting of RI and Exf gw , when using MVR package, is not possible. The RI and Exf gw represent two different physical processes, which should have two different water balance representations-for example, their separate representation can be critical in simulating agricultural irrigation and contaminant transport. After reporting that problem to the US Geological Survey (USGS) team of MODFLOW 6, they confirmed that the RI and Exf gw splitting will be handled in the upcoming version of MODFLOW 6.
In this study, the total streamflow is defined as q = RI s + Exf s gw + q B , where RI s is rejected infiltration routed to streams; Exf s gw is groundwater exfiltration routed to streams; and q B is baseflow. Applying the CRR concept as proposed in this study, the Hortonian flow (q H ) and the Dunnian flow (q D ) are inherently simulated through RI s and Exf s gw components. Fig. 10 shows two cases with two different streamflow states: (1) Hortonian flow (q H = RI s ), also known as infiltration excess runoff; and (2) Dunnian flow (q D = RI s + Exf s gw ), also known as saturation excess runoff. In principle, there is a possibility to define q H and q D spatiotemporally from the MODFLOW 6 output files, but it needs additional script (e.g. Python script) to differentiate spatially and temporally between the areas that will have q H and the areas that will have q D . As such differentiation does not affect water balances, it was not implemented in this study.

Conclusions
The new MODFLOW 6 IHM code was used in transient simulation of surface-water/groundwater interactions in the hard-rock system (HRS) of Sardon catchment (~80 km 2 ) in Spain. The catchment is characterized by a weathered fractured granitic rock system with a shallow water table, dense drainage network, low storage, moderate permeability, large heterogeneity, low human impact, and availability of an automated monitoring network. As that catchment is located in the water-limited Mediterranean environment, its rainfall has a high temporal variability with distinct wet and dry seasons.
1. The new numerical MODFLOW 6 IHM code was applied. The proposed model solution has undergone the following main improvements as compared to previous modelling attempts in that catchment: (1) a new conceptual model, proposed by Francés et al. (2014), was implemented, which is more appropriate for HRS; (2) the P, PET and E I driving forces were improved; (3) a novel CRR concept was implemented in the standard MODFLOW 6 environment; and (4) an irregular Voronoi grid was used, resembling more realistically the complexities of the study area, such as curvatures of the river network, faults boundaries etc., than a standard regular grid. 2. The introduced CRR concept allows one to cascade the rejected infiltration and/or the exfiltrated water between irregular cells, and at each cell, to partition that water between: (1) evaporation; (2) reinfiltration; and (3) stream discharge. The flow partitioning factor (β i,j ), used in the CRR, shows a significant impact on that water partitioning between ET and q and consequently affects the water balance; therefore, control of β i,j , using additional observations, is recommended. 3. The introduction of the CRR concept implied that nearly all WB equations had to be reshaped; the remaining difficulty in that respect was that the RI and Exf gw , when using the MVR package, could not be separated in the current version of MODLFOW 6; once this is overcome, the WB equations, proposed in this study, will become more straightforward. 4. The main outcomes of this study were similar to the outcomes of the earlier GSFLOW study by Hassan et al. (2014), i.e.: (1) large R g and Exf gw due to the shallow water table and moderate permeability, low storage and dense drainage network-enhancing groundwater discharge to streams; and (2) large spatio-temporal variability of groundwater fluxes implying that spatial recharge/discharge conditions conformed to the typical HRS-WLE mosaic pattern of R n . However, as a result of changes listed in section 'Methods', the following differences were notified: (1) the E s substantially increased; (2) the R g and Exf gw increased but proportionally so that the R g -Exf gw difference remained nearly the same; and (3) the R n became more 'modest', i.e. was less positive in the wettest years and less negative in the driest years. 5. With the current capabilities of MODFLOW 6, including the CRR proposed in this study, the single-environment MODFLOW 6 has capacity that is already comparable to the multienvironment IHMs but has the important advantages of larger flexibility and faster processing time.
The following further improvements of MODFLOW 6, bringing it even closer to the complex, multienvironment IHMs, are recommended: (1) implementation of the E s in the MODFLOW 6 environment; (2) improvement of the surface flow component of MODFLOW 6 (direct runoff), e.g. as proposed in this study (CRR), (3) enhancement of the streamflow simulation via the capabilities of the SWR1, or similar; and (4) improvement of the ET g simulation, or, even better, full partitioning and sourcing of ET ss , e.g. as proposed by Francés (2015).