Regional coupled surface–subsurface hydrological model fitting based on a spatially distributed minimalist reduction of frequency domain discharge data

. Although integrated water resource models are indispensable tools for water management at various scales, it is of primary importance to ensure their proper ﬁtting on hydrological variables, avoiding ﬂaws related to equiﬁnality. An innovative stepwise ﬁtting methodology is therefore proposed, which can be applied for any river basin model, from catchment to continental scale as far as hydrological models or land surface models are concerned. The methodology focuses on hydrosystems considering both surface water and groundwater, as well as internal water ﬂuxes such as river baseﬂow. It is based on the thorough analysis of hydrological signal transformation by various components of a coupled surface–subsurface hydrosystem in a nested approach that considers the conditionality of parameter ﬁelds on their input forcing ﬂuxes. The methodology is based on the decomposition of hydrological signal in the frequency domain with the HYMIT (HY-drological MInimalist Transfer function) method (Schuite et al., 2019). Parameters derived from HYMIT are used to ﬁt


25
Given the current climate and anthropogenic trends, water management has become one of the greatest challenges of the 21 st century. Today, less than 1% of the global water stock is readily available for human activities, among which only a fraction constitutes a renewable water stock (Oki and Kanae, 2006;Roche and Zimmer, 2006;de Marsily, 2008). Overall, the pressure on the earth hydrosystem is such that approximately 500 M people are already experiencing water stress throughout the year (Mekonnen and Hoekstra, 2016). Although uncertainties remain regarding the quantification of climate change impacts on 30 water resources (Taylor et al., 2013), the number of people exposed to hydrological stresses will rise in the near future. The main source of uncertainties for such predictions is related to climate models, but uncertainties related to hydrological models are not negligible (Hattermann et al., 2018;Her et al., 2019;Ashraf Vaghefi et al., 2019). It is therefore important to ensure the proper adjustment of these models to hydrological data.
In the Anthropocene Epoch (Crutzen, 2002;Crutzen and Steffen, 2003), it is not only climate that affects the trajectory of 35 large river basins but also human pressure, such as land cover changes and water withdrawal that can sometimes become the controlling factors of those systems, overcoming natural factors (Rockström et al., 2009). In some regions of the world, water demand for sociological purposes will also drive water cycle changes in the same order of magnitude as climate change could affect the system (Haddeland et al., 2014). It is therefore now acknowledged that integrated watershed modeling tools are the most suited for water management and planning purposes (Perkins and Sophocleous, 1999;Flipo et al., 2014;Barthel and 40 Banzhaf, 2016), usually based on scenario testing (Hattermann et al., 2017), as well as for disentangling climate change impacts on the water cycle from a growing water demand associated with population growth (Jackson et al., 2001;Flipo et al., 2021a).
In this context, improvements in hydrological modeling tools are needed in order to properly estimate the spatio-temporal dynamics of water fluxes involved in the terrestrial water cycle (Uniyal et al., 2015). Regarding climate impact studies, it is of utmost importance to ensure the proper calibration of hydrological models before using them as prospective tools. 45 Since the hydrosystem modeling blueprint proposed by Freeze and Harlan (1969); Simmons et al. (2020), and efficiently implemented in many well-known models, hydrology and hydrogeology communities alike have followed divergent paths leading to a paradoxical situation where even though each discipline depends on the other through boundary conditions, collaborative work between both sides needs to be reinforced (Staudinger et al., 2019). As a concrete consequence of this situation, modelers are torn apart between more hydrological interests, usually at the continental scale where the calibration of the models is 50 mostly focused on discharge, and catchment or regional scale where occasionally hydrogeologists add hydraulic head to the fitting process through reformulation of objective functions (Saleh et al., 2011;Flipo et al., 2012;Pryet et al., 2015;Baratelli et al., 2016). All these approaches suffer from large uncertainties in the identification of parameters known as "equifinality" (Beven and Binley, 1992;Beven, 2006;Ebel and Loague, 2006). On the one hand, fitting model parameters, and especially groundwater (GW) models, on discharge data does not prove neither that the model reproduces the correct physical processes, 55 nor that the distribution of river-aquifer exchanges is correctly located along a river network at the watershed scale (Barclay et al., 2020). On the other hand, the most recent benchmarking strategies rely on the ability of the model to reproduce physical processes based on simplistic case studies but not on a data-model comparison strategies (Maxwell et al., 2014;Tijerina et al., 2021) that overcome equifinality issues.
Three main types of fully integrated model structures can be distinguished: (i) fully physically-based 3-D models (e.g. Hy-60 drogeosphere (Therrien et al., 2010), ParFlow (Kollet and Zlotnik, 2003), InHM (VanderKwaak and Loague, 2001;Loague et al., 2005), CatHy (Camporese et al., 2010)), (ii) fully physically based pseudo 3-D models (e.g., Mike-she (Abbott et al., 1986a, b) or (ii) coupled conceptual-physically-based models (e.g., CAWAQS-like (Flipo et al., 2021a), GSFLOW (Markstrom et al., 2008)), for which surface processes are simulated using conceptual reservoir models. They all compute the hydrological processes controlling hydrosystem 1 functioning. They are thus particularly suited for reporting the spatial and temporal dy-65 namics of water fluxes for water management purposes . However, results provided by such integrated models can be uncertain (Wu et al., 2014). These issues can be due to the significant number of calibration parameters involved (Wu et al., 2014) and to the reliance of subsurface parameters on recharge rates estimated by simulation of surface processes (Erdal and Cirpka, 2016). In order to reduce integrated model uncertainties, their parameters need to be defined more precisely, through specific calibration procedures. In this sense, Flipo et al. (2012) introduced a step-wise calibration strategy of hy-70 drosystem models in which surface and subsurface parameters of hydrosystem models are calibrated in a sequential fashion to address their dependency. In this procedure, surface and subsurface models are iteratively optimized until the calibrated parameter set reproduces both observed groundwater levels and river discharges. This procedure introduces the fluxes occurring at the surface-subsurface interface (aquifer recharge and stream-aquifer interactions) in the calibration procedure indirectly and accounts for the dependence of subsurface parameters on surface recharge. However, even if the step-wise calibration strategy 75 has proven its efficiency in fully coupled model calibration (Flipo et al., 2012;Pryet et al., 2015;Baratelli et al., 2016), some aspects remain critical such as the computational burden of conducting the iterative procedure and the potential bias in the simulation of the water budget.
In this paper, we demonstrate a step-wise methodology that builds on Flipo et al. (2012). It considers the conditionality of GW parameters on their boundary conditions and relies on estimates of internal hydrosystem flux values to improve the calibra-80 tion of spatialized hydrosystem parameters. Besides bringing together the interests of the two communities of hydrologists and hydrogeologists, the methodology intends to significantly reduce the equifinality issue related to the fitting of simulation models of hydrosystems, in terms of the watershed. The methodology can be applied for any river basin model from the catchment scale to the continental scale as far as Land Surface Models (LSM, Pitman (2003)) are concerned, such as VIC (Liang et al., 1994(Liang et al., , 1996Hamman et al., 2018), SURFEX (Masson et al., 2013;Le Moigne et al., 2020), ORCHIDEE (Ducharne et al., 85 2003;Krinner et al., 2005), or MGB-IPH (Collischonn et al., 2007;Paiva et al., 2013), or mHM (Samaniego et al., 2010). The methodology focuses on hydrosystems considering both surface water and GW, as well as internal water fluxes such as river baseflow. It is based on the thorough analysis of hydrological signal transformation by various components of a hydrosystem in a nested approach, which takes into account the conditionality of parameter fields on their input forcing fluxes.
The methodology, described in detail in section 2, is based on the decomposition of a hydrological signal in the Fourier 90 frequency domain with the HYMIT (HYdrological MInimalist Transfer function) method (Schuite et al., 2019;Schuite, 2022).
HYMIT makes it possible to study the influence of the physical properties of watersheds on the deformation of precipitation signals. It is based on a transformation of hydro-meteorological data in the Fourier frequency domain, in which a transfer function is composed of a minimalist number of hydrological processes and parameters are fitted with a Markov chain Monte Carlo (MCMC) approach. Parameters derived from HYMIT are used to fit the coupled surface-subsurface hydrological model 95 CaWaQS3.02 with a step-wise approach, which relies on successive MCMC optimizations, respectively, related to various objective functions that represent the dependency of hydrological parameter fields on forcing input fluxes.
While section 3 explores the performance of the model from both quantitative and qualitative perspectives, section 4 illustrates the significant progress in terms of spatialization of the water balance components, as well as stream-aquifer exchanges, at the regional scale, enabled by the proposed methodology. It is exemplified with the Seine River basin (∼ 76,000 km 2 ).

100
More specifically, these two sections reveal the consistency of the spatially distributed physical properties of the Seine subcatchments, derived from the reproduction of the partitioning of effective rainfall into fast surface flow and slow flow in the subsurface domain, with morphological data and geomatic analyses. They highlight the evolution of the filtering effect of precipitation signals by successive catchments from upstream to downstream at the scale of the Seine basin. Finally, a discussion (section 5) points out the relevance of the developed methodology for hydro(geo)logical modeling in general, from the 105 catchment to the continental scale.
2 Material and method: Fitting a regional river basin hydrosystem model Fitting a regional hydrosystem model, which couples fast surface and slow subsurface flows at the scale of a regional basin, such as the Seine River basin, is a challenging issue (Flipo et al., 2012), especially regarding GW, since measurements of the fluxes are unavailable.

110
Although automatic adjustment methods do exist in disciplinary fields, whether in hydrology or hydrogeology, this is not the case when one is interested in the coupled dynamics of surface and underground processes. In their review of this topic, Flipo et al. (2012) proposed to circumvent the problem using a nested loop fitting method. The basis of the method is to accept the conditionality of GW flows on the process of aquifer recharge through its impluvium. It is then possible to nest two loops: one dealing with fast surface processes and the other associated with slow subsurface processes. An essential condition for the implementation of this method is that the first loop enables the quantification of the forcing of the second loop, i.e., recharge of aquifers by their impluvium.
An initial implementation of this two-step methodology, each one being fully automated, was proposed by Labarthe (2016), and has been used to reconstruct the hydrological trajectory of the Seine basin since the beginning of the XX th century (Flipo et al., 2021a). This method, relying on the fact that, for a 17-year stationarity period, river baseflows can be equated 120 to the aquifer recharge across the river basin, depends on additional information to the classic datasets of river discharges and groundwater hydraulic heads, namely, an estimate of the river baseflow. However, these estimates are still marred by an 4 https://doi.org/10.5194/gmd-2022-24 Preprint. Discussion started: 14 April 2022 c Author(s) 2022. CC BY 4.0 License. error that we acknowledge, but that we are unable to quantify, mainly due to the great difficulty (or even impossibility) of the measurement at large scale.
We therefore develop a new fitting methodology based on a thorough analysis of hydrological signal transformation by 125 various components of a hydrosystem in a nested approach, which considers the conditionality of parameter fields on their input forcing fluxes. The approach is based on the decomposition of hydrological signals in the Fourier frequency domain with the HYMIT method (Schuite et al., 2019). The Seine River basin ( 76,000 km 2 ) is a highly human-impacted basin with a population of 17 M inhabitants. It receives the highest anthropogenic pressure in France, due to the industry and agriculture linked to the development of the urban area of Paris, which coexists today with highly productive agricultural areas (66 % of the basin area, Fig 1). For more information on the Seine basin, please refer to Billen et al. (2007) and Flipo et al. (2021b).
The river network is composed of 28,000 km of perennial rivers (Fig. 2). Overall, 97 % of the river network lies within the 135 sedimentary Paris basin (Guillocheau et al., 2000), the largest GW reservoir in Europe. Water withdrawals in surface water and GW amount to 3 km 3 a −1 . The pressure on water resources is high especially during 140 low-flow periods, when river discharges are sustained by GW and by human-built reservoirs (Fig. 2), totalling 841 10 6 m 3 .

Naturalization of downstream discharges from reservoirs
Here, observed river discharges are studied as an overall integrator of the system's behavior and response to climatic fluctuations. In this context, preliminary corrections are required to ensure the absence of anthropogenic disturbances in measured data that might result from the significant pressure the basin is submitted to. Since the HYMIT approach is not designed to 145 dissociate such perturbations in measured discharge variations, the imprints left in the observed data by water reservoir storage/release cycles have been subtracted from downstream station records. Two types of information were used to perform this discharge naturalization: (i) reservoir withdrawn or released volumes to the river system and (ii) water travel times along the network, respectively, associated with each downstream station from one or several reservoirs.
The calculation of a naturalized discharge value Q i nat (t) [m 3 s −1 ] associated with its respective actual measurement Q i mes (t)

150
[m 3 s −1 ] is performed using equation (1), written in the case of a station i, at a time t, located downstream from N reservoirs: where Q i dam,k (t − T i tra,k ) [m 3 s −1 ] represents the daily volume either stored from (> 0) or released to (< 0) the river system, accounting for the water travel time T i tra,k [s] along the network fraction between dam k and station i. Travel times are calculated for each reach r of the network, using a relative transfer time index I tr (r) = dl(r)/ s(r) S(r) γ , 155 as a function of geomorphological data (Golaz-Cavazzi, 1999;Flipo et al., 2012): the distance dl(r) [m] between center of reach r and its contiguous downstream reach; reach slope s(r) [−], derived from DEM; the cumulative reach upstream drainage area S(r) [km 2 ] and a calibration parameter γ = 0.25 (Korkmaz, 2007).
For every reach, a relative transfer time index to the basin outlet I tr;r→out is calculated as a sum of local I tr along all reaches 160 leading to the outlet. Expression (2) is then used to characterize travel time T i tra,k , where I tr;max is the maximum relative time  index to the outlet basin and T c [s] is the basin global concentration time, considered equal to 17 days (Gomez, 2002;Saleh et al., 2011;Pryet et al., 2015).
Naturalization of discharge records was carried out for the 30 measurement sites mentioned in figure 2, for which travel time 165 values from reservoirs to stations (in days) are also displayed. In the case of several upstream reservoirs, labels correspond to mean travel time values.  conductance coefficient, Qenr: stream-aquifer exchanges flow, Q lim : river-to-aquifer limit infiltration rate, ∆ h: difference between water height in river (hriv) and aquifer hydraulic head (haq).
The physically-based CaWaQS coupled model (CAtchment WAter Quality Simulator - Flipo et al. (2005Flipo et al. ( , 2007bFlipo et al. ( , a, 2021a; Labarthe (2016)) was used to: (i) estimate the distributed effective rainfall over the basin as an input of a distributed HYMIT  (1978) and implemented as the MODCOU-NEWSAM software suite (Ledoux et al., 1984(Ledoux et al., , 1989, CaWaQS3.02  is a spatially distributed model that simulates coupled water, matter, and energy balances and flow dynamics within all compartments of a hydrosystem. The software structure links dedicated C-ANSI libraries meant to mimic main physical processes controlling the fate of water in 175 each compartment. Calculations of surface, subsurface, and GW dynamics involve five main modules using a daily time-step (Labarthe, 2016) : a surface module (Fig. 3a), which computes estimates of AET, runoff, and infiltration fluxes. Relying on a conceptual reservoir-based approach (Girard et al., 1980;Deschesnes et al., 1985), water balance calculations account for climate data (see rainfall and PET maps in Figs. 8a and 8b, respectively) as well as the distributions of land use and parent soil 180 material ( Fig. 1). Runoff water production is aggregated according to local subcatchments to ensure its direct transfer to the river system. These catchments can also be set up as runoff short-circuits toward the unsaturated zone ("chasm"-type configuration); a vadose zone module (Fig. 3c), which vertically transfers infiltration from the surface domain to subsurface outcropping aquifer areas. It diffuses soil infiltration based on a Nash reservoir cascade (Nash, 1959;Besbes and De Marsily, 1984) 185 toward the aquifer system; a groundwater or aquifer system module ( Fig. 3d), based on the pseudo-3-D diffusivity equation (de Marsily, 1986), solved using a semi-implicit finite volume numerical scheme, applied on nested grids. Besides integrating water recharge and anthropogenic withdrawals, it accounts for both confined-and unconfined-related resolution particularities, and also handles, along time and space, reversible transitions between these two states. Exchanges between aquifer units are 190 simulated based on a 1-D vertical simplification of water fluxes, which are assumed to be linearly connected to the head difference between aquifer units; a non-linear conductance model (Fig. 3e), which accounts for a limitation of the infiltration flux in the case of disconnection (Brunner et al., 2009;Rivière et al., 2014;Newcomer et al., 2016), is integrated within a Picard-iterative approach to compute stream-aquifer exchanges (Rushton, 2007;Ebel et al., 2009;Flipo et al., 2014);

195
a hydraulic module (Fig. 3b), which transfers in-stream water discharges using a Muskingum routing scheme (David et al., 2011(David et al., , 2013. For each river network cell, computed discharges integrate stream-aquifer fluxes, inputs due to subsurface runoff as well as exogenous point injection flows. a river network that corresponds to 6,830 km of rivers. Mainly due to computational time concerns, calculations of stream-aquifer exchanges have been constrained to rivers including reaches with Strahler orders (Strahler, 1957)  The method requires knowledge of the effective rainfall, which is not directly measured. It therefore needs to be assessed by a surface balance model. In this study, the surface module of the hydrological model CaWaQS3.02 is used (see section 2.2) 240 both for practical reasons (spatialized soil, land surface data availability, and integration) and for consistency concerns with subsequent fitting steps. The fitting of the surface balance model is detailed in the next section 2.4.
The power spectral density (PSD) of the effective rainfall sub-basin average is taken as input for the HYMIT parameter estimation. The PSD of the naturalized river discharges at the outlet is calculated as the output. The ratio of the latter over the former gives an experimental transfer function (ETF) (Pedretti et al., 2016;Schuite et al., 2019). An MCMC inversion 245 procedure is implemented to adjust HYMIT to each ETF in order to estimate all five parameters for all sub-basins for which discharge data are available. The MATLAB package of Haario et al. (2006) for MCMC inversion is used with the built-in Metropolis-Hastings sampler.
A total of 384 gauging stations are (or were) operational across the different rivers of the Seine basin. Daily discharge data were checked for completeness and overall quality. In the presence of any gap longer than 10% of the total series length, 250 the station was discarded from the analysis. So were stations with overly short records (less than 10 years), to maximize the statistical power of the analysis. Small gaps are filled by linear interpolation. Minor data errors were corrected, such as inconsistent timelines, double records, or negative values. After data curation and selection, discharge series from 221 stations remained. Prior to Fourier transformation, all time series are detrended and windowed using a Hanning tapering function.

2.4
Step-wise fitting methodology based on a nested hydrological approach 255 A nested fitting method is proposed (Fig. 5). It is based on the identification of the CaWaQS3.02 model parameters, process by process, on the basis of measurable quantities. It is considered nested due to the conditionality of the parameter fields on their forcings. Each step is developed with this conditionality idea in mind. Roughly, the first step considers the estimation of the partitioning of the rainfall into AET and effective rainfall, while the next one focuses on the partitioning of effective rainfall into fast runoff and slow infiltration. The last steps deal with regulating the velocities of both surface and subsurface flows.

260
The first step of the fitting methodology consists in estimating the total AET at the basin scale from the average discharge of water flowing out of the basin (Fig. 5, step 1). In order for this average to exist mathematically, it is then necessary to reproduce the flows averaged over 17 years, which is the stationarity period of this signal (Flipo et al., 2012;Massei et al., 2010Massei et al., , 2017. In the surface balance module, the set of parameters regulating actual evapotranspiration (AET) flow simulation are adjusted using an MCMC approach aiming at minimizing the discrepancy between the long-term average discharge at the 265 watershed outlets and the long-term average meteoritic net input over the watershed impluvium (i.e., the effective rainfall that corresponds to water available for flow within the basin). Furthermore, we selected a total of 35 reference sub-basins across the Seine hydrosystem satisfying two criteria, namely, good spatial coverage of the basin as well as long and complete daily discharge time series.
Before proceeding with the fitting of the model parameters themselves, sets of distributed minimalist hydrological parameters 270 are estimated at 221 gauging stations in the basin using the HYMIT method (Fig. 5, step 2). The following steps of the fitting procedure are based on these estimates with the objective each time to spatially reproduce the variability of these minimalist parameters with the process-based model used, e.g., CaWaQS. In other words, the wealth of information provided by the quantified and regionalized estimates of HYMIT parameters is used as the only support for the fitting of the surface and subsurface water flow calculation within the model. Thus, this step enables the transition from a point assessment to a 275 continuous spatio-temporal characterization of the regional hydrosystem behavior. The third step consists in adjusting a parameter of the water production cells in order to reproduce as closely as possible the effective rainfall partitioning evaluated by HYMIT (β coefficient, Fig. 5, step 3).
The fourth step deals with the temporal synchronization of fast runoff estimated by the process-based model using the estimate of this quantity obtained with HYMIT at each gauging station (λ coefficient, Fig. 5, step 4). The attribution of λ coefficient 280 values identified with HYMIT to the corresponding CaWaQS parameter is made by analogy. For each gauging station, a spatial analysis is performed to calculate the proportion of the CaWaQS hydrological units that comprise the catchment of the station.
Over the Seine basin, it is possible to identify a sufficient number of upstream catchments that are mostly composed of a single CaWaQS hydrological unit. For each of these catchments, the CaWaQS parameter is set up in a straightforward way, simply equating it with the value of the HYMIT λ coefficient.

285
Until now, all the steps were focused on surface processes. The proper simulation of these processes is crucial since they allow us, among other processes, to estimate a distributed recharge toward the aquifer system. The fifth and last step of the methodology is carried out differently for the vadose zone than for the aquifer system and is divided into three substeps : a preliminary coarse calibration of transmissivity parameter fields for each model aquifer layer, performed in steady state. A trial-and-error approach is used to adjust simulated water tables to mean level values measured at control wells.

290
To achieve this step, an aquifer compartment-only simulation is implemented. A mean aquifer recharge field, calculated over two 17-year cycles (i.e., 1986-2020 period), derived from the HYMIT-calibrated CaWaQS surface module, is used to constrain the simulation (Q i on Fig. 5, step 5); a fitting of the Nash-cascade parameters, namely, k and N (see sections 2.3 and 2.2). Usually associated with lithology, HYMIT-k values are spatially distributed along a functional sectorization combining two information types: dominant 295 soil textures (Fig. 1) and upmost free aquifer formation lithology (Fig. 3). For each association identified, attribution of k values is made using an analogy method, as previously described in the case of HYMIT-λ coefficients distribution (see section 2.4). N parameter field, reflecting unsaturated zone thickness, is geometrically defined, based on the difference between subsurface adjusted water levels in steady state and ground elevation, while considering a single reservoir to be representative of an elementary 5-m thickness. Overall, at the basin scale, k parameter values range from 2 to 9 days.

300
Tertiary terrains are mostly associated with low thickness values, from 0 to 20 m, while ranging up to 40 m within the Chalk impluvium limits (Fig. 11). An unsaturated thickness map ( The next two sections exemplify the power of the step-wise methodology applied to the regional scale of the Seine basin.

Qualitative analysis of spatially distributed infiltration fluxes estimated with HYMIT
The average partitioning of effective rainfall between surface runoff and deep infiltration is adjusted based on the HYMIT β parameter. This parameter is generally poorly constrained at large scale, especially its spatial distribution. Hence, very little information is available to compare our β estimates against, and ultimately, to validate them.

330
Fortunately, the French Geological Survey (BGRM) has developed a systematic method to qualify the propensity of terrains to either infiltrate water or to generate runoff, as described by Mardhel et al. (2021). The method takes advantage of the mismatch evaluated between thalweg locations, inferred from digital elevation models, and the actual location of drainage networks, used to build a normalized index called IDPR (Network Development and Persistence Index). The lower its value, the more a terrain is prone to infiltration and vice versa. Mardhel et al. (2021) calculated this index for the entire French 335 metropolitan area at high resolution (25 m). We take this opportunity to qualitatively assess the consistency of β estimates, which, in a way, is very complementary to the IDPR index, being much less spatially resolved but providing an actual operable value to the flow partitioning, which is of primary importance for modeling purposes. First, an overall spatial coherence is observed between these two indicators. Upstream stations experience medium-to-high where low values of β are found. We further note a satisfactory sectorial coherence in β estimates: Two low-order rivers draining the same geomorphological units exhibit close properties in terms of flow partitioning. This aspect is particularly well established in the "Beauce" and "Craie Champenoise" regions ( Fig. 6, panels (b) and (c), respectively). It is also the case for 345 the southern regions of the Morvan (see Fig. 2); yet in this case β values seem comparatively low in regard to the high IDPR values, possibly underscoring a more complex interplay between different runoff-generating mechanisms in this sector (slope, soil structure, geology, etc.). with high slopes, extending NW-SE (see Fig. 2), among Upper Cretaceous carbonate subplanar units (see Fig. 4). Therefore, property contrasts between these terrains are clearly distinguishable in IDPR values and are also well captured by the flow partitioning estimation with the HYMIT, as the river draining the "Pays de Bray" anticline displays a β value of approximately 45 % to be compared with the very low values found on adjacent stream networks (β < 10 %).

Simulation of river discharges
Similarly to raw piezometry data, a pre-processing work is carried out on river discharge observations using data compiled 375 from the HYDRO database (http://hydro.eaufrance.fr/), totalling 384 river stations (see section 2.3). As previously stated, since calculations of stream-aquifer exchanges are constrained to the main river network (see section 2.2 and Fig. 4), 167 gauging stations only are considered as valid discharge calibration control points. Usual Nash-Sutcliffe (Nash and Sutcliffe, 1970) and respectively,in Tab. 2) are selected to evaluate simulation performances on river discharges.

380
At the basin scale, 50 and 53 % of control stations, respectively, show N SE and KGE scores above or equal to 0.5. The spatial distributions of both criteria (not shown) demonstrate that most of the higher values are noticeably distributed along the Seine River and its eight main tributaries (see bold blue lines in Fig. 2). For the 64 stations located along this network portion (2,910 km, 43 % of total modeled network), 54 (resp. 26) of them are associated with KGE ≥ 0.5 (resp. 0.7). Noticeably, 44 stations combine both N SE and KGE values ≥ 0.5.
The effective rainfall is highly contrasted across the Seine basin, ranging from 44 mm a −1 up to 918 mm a −1 locally 395 (Fig. 8d). The lowest effective rainfall rates occur around the city of Paris and in the Eure basin, around Chartres and Dreux (see Fig. 2). The central part of the basin experiences an effective rainfall mainly lower than 130 mm a −1 , while the eastern ridge, or Jurassic edge, experiences effective rainfall rates higher than 250 mm a −1 , sometimes reaching 920 mm a −1 as is the case in the southern part of the basin in the Morvan area. It is also the case in Normandy, in the northern part of the basin, called "Pays de Caux" (see Fig. 6, panel (a)), in the north of Rouen.

Distribution of infiltration rates
A key piece of information for hydrogeologists and groundwater managers is the estimation of aquifer recharge. A complete water balance calculation using the CaWaQS model enables the simulation of effective rainfall, runoff, and infiltrated water distributions at the daily time-step with a high spatial resolution. It is therefore possible to represent the average distributed partitioning of effective rainfall at the basin scale (Fig. 9), over the period 1970-2018, expressed as the ratio between the 405 simulated runoff and the total effective rainfall. Before analyzing further the estimated infiltration rates, it is important to note that the distributed effective rainfall partitioning estimated by CaWaQS-Seine is consistent with that obtained using the HYMIT analysis alone (Fig. 9), which is in that agree with HYMIT low partitioning values for head water streams in the area (lower than 0.3, Fig. 9).
Distributed infiltration rates are then calculated, as an inter-annual average over the simulation period 1970-2018 (Fig. 10).
They balance the effective rainfall mostly leading to relatively moderate infiltration rates (lower than 120 mm a −1 ) over the basin, except (i) on the Normand "Pays de Caux" north of Rouen where it can reach 400 mm a −1 , and (ii) to a lesser extent  There is also a spatial coherence between β and the infiltration rate from the point of view of land use, insofar as the predominantly agricultural (resp. urban) areas are indeed marked by higher (resp. lower) infiltration (see, for instance, the 420 Chalk area or the central Île-de-France that surrounds Paris city). Also, the areas initially associated with very low IDPR values (fast infiltrating areas -see Fig. 1) clearly appear on the infiltration map, with infiltration rates higher than the regional average (see, for instance, the Pays de Caux area or the eastern of Chartres).
Finally, geological characteristics also emerge from the zoning, insofar as the sectors dominated by soils classically more permeable and conductive to infiltration (e.g., sands, alluvium) delimit areas of higher infiltration rates than the regional trend.

425
Such configuration is noticeable for regions such as the Sologne or Fontainebleau forest areas (see Fig. 1), or Thanetian sands and alluvium deposits (see Fig. 4).

Distribution of groundwater contribution to river discharges
Labarthe et al. (2015) and Pryet et al. (2015) were the first to publish a spatially distributed evaluation of stream-aquifer exchanges at the regional scale, exemplified with the Seine basin. At the time, they followed the step-wise fitting methodology 430 of Flipo et al. (2012). Given the advances proposed in the current paper, we re-assessed those estimates (Fig. 11) based on the HYMIT functional analysis of the hydrological behavior of many sub-basins and also on the extension of the subsurface domain that is taken into account by the model.
The pluri-annual average (2003-2020) contribution of groundwater (GW) to river discharge is calculated along the river network of the Seine basin. Two specific patterns of spatial GW contribution to river baseflow appear in the Seine basin:

435
a longitudinal increase in the contribution of GW to river baseflow from upstream to downstream for river systems originating in the Jurassic edge of the basin; a very high contribution of GW to river baseflow over the tertiary (> 0.75), alluvial (> 0.4), or Upper Cretaceous (mostly > 0.6) aquifer units (Fig. 11, panel (b)).
Those patterns are confirmed by a spatial analysis of the statistical distribution of GW contribution to river baseflow regarding 440 the specific discharge. Each analysis is made for each outcropping aquifer unit. For the Jurassic and Lower Cretaceous aquifer units, the lower the specific discharge, the higher the GW contribution to river baseflow with a highest limit at 0.7 for a few small streams in the area (Fig. 11, panel (b 4 )). For the small streams fed by GW from the Tertiary aquifer units, the GW contribution is always higher than 0.8 (Fig. 11, panel (b 2 )). The GW contribution of the Upper Cretaceous aquifer unit is similar to that of the Tertiary aquifer units (Fig. 11, panel (b 3 )), even though the absolute value of the contribution can be 445 slightly lower in this area (0.5). As for the Jurassic and Lower Cretaceous edge of the basin, the GW contribution of alluvial GW to river discharge exhibits a decreasing GW contribution with the increase of the specific discharge that is correlated to higher Strahler orders (Fig. 11, panel (b 1 )). But the GW contribution in alluvial aquifers is always higher than 0.4, which is a significant difference with the Jurassic edge. Finally, a water budget of the Seine basin is established (Fig. 12) over a 17-year cycle that ensures zero storage variations (Flipo et al., 2012).  Anthropogenic withdrawals from the aquifer system account for 7.8 % of the total aquifer recharge (approx. 1.09 km 3 a −1 ).
As discussed in section 4.3, the exfiltration regime from the aquifer system to the hydraulic network largely dominates river-aquifer exchanges, as shown by the positive net exchange values in figure 12. Within the basin limits, flows drained from the 460 underground system are, on average, responsible for 67 % (318 m 3 s −1 ) of the Seine River discharge at the outlet of the basin (477 m 3 s −1 ). A complementary input to the river network is composed of runoff contributing to 160 m 3 s −1 , which accounts for 33 % of the river discharge at the outlet of the basin.
One main advantage of distributed models is the ability to calculate water budgets also for sub-basins of the regional Seine basin. More details about the spatial distribution of the water budget over the Seine basin are provided in Tab. 3.

465
5 Discussion on the relevance of the proposed methodology across scales The step-wise fitting methodology depends on a distributed minimalist reduction of hydrological parameters, which is performed using the HYMIT methodology (Schuite et al., 2019) coupled with optimization of the distributed parameters of the forward model CaWaQS3.02 using mostly MCMC optimizations or set up by analogy. HYMIT is a very powerful method that may provide a consistent view of minimalist hydrosystem hydrological functioning at various scales given that its fundamental 470 hypotheses are fulfilled; therefore, it is a useful companion for adjusting various hydrological models from the catchment scale to the continental scale (Flipo et al., 2014).

Fundamental hypotheses behind the step-wise methodology
The proposed methodology relies on two main assumptions that have to be fulfilled before it can be applied: the stationarity of hydrological signals over the period of time for which the study is run, especially for performing the 475 transfer function analysis with the HYMIT, which involves effective rainfall and river discharges; the overlapping of surface and subsurface catchments.
Before discussing these hypotheses and their validity across scales, let us introduce the instantaneous equation for the water budget at the surface basin scale: The first step in the step-wise methodology relies on the estimate of aet, which is done by simplifying Eq. (3). The two 485 hypotheses make this step possible on whatever scale it is applied: from small catchment scale to continental scale. The second step of the step-wise methodology also relies on these two assumptions.

Overlapping of surface and subsurface catchments
The inflow or outflow through the limits of the surface topographical basin projected onto the subsurface groundwater domain is a quantity that is neither observed nor measured, whatever the scale of the system of interest, making the estimation of the 490 term q bound (t) of Eq. (3) almost impossible. Assuming that the surface and subsurface watersheds overlap (Tóth, 1962), the subsurface fluxes at the basin boundaries (lateral and bottom) can then be neglected (i.e., q bound (t) = 0).
This hypothesis holds in sedimentary basin environments where fractured or karstified areas are not preponderant. Tóth (1962) shows that in sedimentary basins, most of the river water fluxes originate from shallow subsurface flows and that piezometric heads are strongly correlated with surface topography. Thus, on a sedimentary hydrosystem, all the watersheds of 495 gauged tributaries of the main hydrographic network can be considered as sub-hydrosystems.

Stationarity of hydrological signals
From a theoretical standpoint, signal stationarity is ensured when its mean and variance remain constant over time. The presence of an underlying trend in hydrological records causes a variation in mean, whereas multi-scale natural fluctuations may affect variance stability. Signal stationarity is particularly important for the two first steps of the nested fitting approach presented in 500 this paper.
Large-scale climate oscillations induce large periodic variations of both surface and subsurface water stock and fluxes (Flipo et al., 2012;Massei et al., 2010). In the absence of long-term trends, these quantities are stationary over major pluri-annual hydro-climatic periods, such as the North Atlantic Oscillation (Flipo et al., 2012). Over such a climate period, the time integrals of storage variations in subsurface and surface compartments are therefore negligible ( ∆s sub (t) dt ≈ 0 and ∆s riv (t) dt ≈ 505 0).
Integrated over a major hydro-climatic period, Eq. (3) yields: where P [m 3 s −1 ], AET [m 3 s −1 ] and Q out [m 3 s −1 ] are the averaged precipitation rate, actual evapotranspiration, and total discharge at the basin outlet, respectively.

510
The first step in the methodology hence requires the identification of stationary time windows for any given hydrosystem, given the a priori knowledge of its pluriannual modes of climatic variability (in the case of the Seine basin, 17 years). Over such a period for which hydrological signals are stationary, AET can therefore be estimated from classic hydrological data that are either measured in situ or spaceborne.
For the second step, namely, the application of HYMIT, first-moment stationarity is routinely forced by detrending signals 515 prior to frequency-domain transformation. By taking into account only the longest and most complete hydroclimatic datasets, which is necessary in order to benefit from the statistical power of HYMIT analysis, we also maximize variance stability over time, owing to the large ratio of short-term fluctuations to long-term oscillations.

Estimating hydrosystem inner fluxes across scales
We developed herein a step-wise fitting procedure that leverages the HYMIT method and that displays unprecedented performance especially for simulating hydrosystem inner fluxes, which are the most uncertain in many current calibrated/validated hydrosystem models, either Land Surface or Hydrological Models (LSMs and HMs) .
One explanation could be the equifinality that stems from large uncertainties in the identification of parameter values (Beven and Binley, 1992;Beven, 2006;Ebel and Loague, 2006). Not only the distribution of actual evapotranspiration or soil moisture are potentially miscalculated (Stisen et al., 2011;Rakovec et al., 2016) but also the distribution of river-aquifer exchanges 525 remains uncertain along a river network at the watershed scale (Barclay et al., 2020). It is therefore of primary importance to further develop fitting methodology for HMs and LSMs (O'Neill et al., 2021).
Although it remains important to check the model ability to reproduce physical processes based on a simplistic case study (Maxwell et al., 2014;Tijerina et al., 2021), it is still very preliminary in terms of model development and not sufficient to meet the challenge of the hyper-resolution that imposes hydrological predictions to be relevant "everywhere" on earth, either 530 at the outlet of large river basins or at the local catchment scale of a few ha or km 2 (Wood et al., 2011). Using multi-scale basin outlets over a continental scale as a basis of the objective function improves LSM performance significantly but does not prevent them from providing divergent results for various resolutions . It, nevertheless, narrows the issue of equifinality. Also, introducing supplementary hydrological fluxes in addition to river discharges to objective functions usually improves model performance (Baroni et al., 2019). HYMIT has the advantage of providing multi-scale estimates of 535 hydrosystem parameter values that are the basis of the estimation of inner fluxes such as regional-scale runoff expressed by beta times the effective rainfall. It provides such estimates in a fully physically consistent framework that is based on a Fourier domain minimalist reduction and therefore provides invaluable additional data to shape objective functions that reduce equifinality drastically.
Important findings about equifinality were recently reported. First, Cuntz et al. (2015) argued that equifinality may be an 540 artifact of flawed calibration procedures that focus on non-sensitive parameters and unsuitable objective functions. Just as important are the findings of , who show that most of the state-of-the-art LSMs and HMs do not fulfill flux-matching conditions across scales. As a consequence, they do not have consistent hydrologic parameter fields across scales and more problematic is that their parametrization at large continental scale is still unresolved. One way to overcome this issue is then to use multiscale parameter regionalization (MPR), as proposed by Samaniego et al. (2010) and further improved by 545 Kumar et al. (2013), for LSMs as successfully demonstrated by Mizukami et al. (2017) and . MPR relies on the identification of intrinsec hydrogeophysical parameters on which scaling operators are based to ensure flux continuity over scales. Combining this powerful calibration technique with the estimates of inner fluxes of hydrosystems provided by HYMIT could be the next step toward improving the robustness of LSMs and HMs and therefore their prospective power for assessments of climate change impact. The nested hydrological fitting step-wise methodology that assumes the dependency of 550 parameter fields on boundary conditions, especially boundary conditions of fluxes, that was used on the Seine basin could be adapted to identify intrinsic parameter fields such as formulated in mHM rather than HRU-based parameters as is the case in CaWaQS3.02 for the surface compartment.
Finally the value of the step-wise methodology was demonstrated at the regional scale. As discussed above, its potential at the continental scale, which can be viewed as a collection of regional systems (Flipo et al., 2014), seems important and needs to 555 be tested. It would also be challenging to evaluate in a heuristic way until which fine scale the hypothesis of the overlapping of surface and subsurface catchments holds. The results of this evaluation could provide a breakthrough in hydrological modeling of the critical zone, since most of the data are acquired in fine-scale catchments, for instance, at the scale of critical zone observatories (Gaillardet et al., 2018). The Orgeval catchment located in the Seine basin corresponds to a hydrosystem on which (i) many long-term high-frequency hydrological datasets exist (Mouhri et al., 2013;Floury et al., 2017) and (ii) such 560 an evaluation could be carried out especially to elucidate water pathways that constitute fundamental information for the understanding of the biogeochemical behavior of hydrosystems (Floury et al., 2019;Tunqui Neira et al., 2020).  (Flipo et al., 2021a) at both the regional and the territory scales.

Conclusions
The methodology is exemplified with the Seine basin, offering for the first time a very detailed picture of the basin functioning as a whole but also at the scale of territories. All hydrological fluxes are estimated in a consistent way between the frequency and the time domains, from classic ones such as AET, river discharges, to more challenging ones such as GW con-570 tribution to river discharge, or exchanges between aquifer units, via the share between fast runoff and slower infiltration. To the authors' knowledge, it is the first time that such a characterization of hydrological system behavior of a regional-scale river basin is proposed, and which reproduces the observations fairly well in both time and frequency domains.
This development paves the way for significant breakthroughs in hydrological modeling of systems over a large range of scales, from small catchments to regional/continental river basins.

575
Code and data availability. CaWaQS3.02 is available under an Eclipse Public Licence 2.0 in the following zenodo deposit . MATLAB-based HYMIT scripts are available under an Eclipse Public Licence 2.0 in the following zenodo deposit (Schuite, 2022).