# Evaluating Richards Equation and Infiltration Capacity Approaches in Mesoscale Hydrologic Modeling

A. Kholis, T. Kalbacher, O. Rakovec, F. Boeing, M. Cuntz, L. Samaniego

► To cite this version:

A. Kholis, T. Kalbacher, O. Rakovec, F. Boeing, M. Cuntz, et al. Evaluating Richards Equation and Infiltration Capacity Approaches in Mesoscale Hydrologic Modeling. Water Resources Research, 2025, 61 (8), pp.e2024WR039625. 10.1029/2024WR039625. hal-05265611

HAL Id: hal-05265611

https://hal.inrae.fr/hal-05265611v1

Submitted on 17 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2024WR039625

**Key Points:**

* The multiscale parameter regionalization technique enables the derivation of transferable parameters for the one-dimensional Richards equation across varying scales and locations

* The hydrologic model mHM with RE (mHM-RE) preserves robust streamflow performance and enhances the physical consistency of soil moisture (SM) dynamics, albeit with increased computational demands

* The mHM model with the infiltration capacity approach (mHM-IC) provides a computationally efficient alternative for streamflow simulation and reproducing average SM anomalies, though it is less accurate in deeper soil layers

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Kholis and L. Samaniego, afid.nur-kholis@ufz.de; luis.samaniego@ufz.de

**Citation:**

Kholis, A., Kalbacher, T., Rakovec, O., Boeing, F., Cuntz, M., & Samaniego, L. (2025). Evaluating Richards equation and infiltration capacity approaches in mesoscale hydrologic modeling. *Water Resources Research, 61*, e2024WR039625. https://doi.org/10.1029/2024WR039625

Received 11 DEC 2024
Accepted 12 JUL 2025

**Author Contributions:**

**Conceptualization:** A. Kholis,

T. Kalbacher, L. Samaniego

**Data curation:** A. Kholis, T. Kalbacher, O. Rakovec, F. Boeing, M. Cuntz

**Formal analysis:** A. Kholis

**Funding acquisition:** T. Kalbacher

**Investigation:** A. Kholis

© 2025. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Evaluating Richards Equation and Infiltration Capacity Approaches in Mesoscale Hydrologic Modeling

A. Kholis<fnref n="1" />,<fnref n="2" /> ID, T. Kalbacher<fnref n="1" /> ID, O. Rakovec³,⁴ ID, F. Boeing<fnref n="2" />,³ ID, M. Cuntz⁵ ID, and L. Samaniego<fnref n="2" />,³ ID

**Abstract** This study compares two widely used approaches for modeling soil moisture (SM) infiltration in mesoscale hydrology: the one-dimensional Richards equation (1D RE), which governs vertical flux exchange but is nonlinear and computationally demanding, and the infiltration capacity (IC) scheme, which is simpler and restricts SM movement to the downward direction. A major challenge in implementing the RE is the estimation of effective parameters at the typical model resolution (hundreds to thousands of meters), as the equation was originally developed for finer scales. To address this, we conducted experiments using the mHM model with Multiscale Parameter Regionalization (MPR) to parameterize both RE and IC approaches. The RE parameterization relied on three distinct pedo-transfer functions (PTFs). Parameters were estimated across 201 basins in Germany and evaluated using streamflow data at multiple spatial resolutions, as well as in situ SM observations from 46 sites (0–25 cm) and 42 sites (25–60 cm and 0–60 cm). Results show that both mHM-IC and all mHM-RE variants perform comparably in streamflow prediction. The use of MPR enables the spatial transferability of PTF-based parameters. Owing to its two-way flux capability, the mHM-RE variant better captures SM dynamics, particularly in deeper soil layers. Although the IC scheme often leads to saturation at depth, it still provides strong performance in capturing SM anomalies. Overall, the study demonstrates that with appropriate parameterization, the RE approach can yield transferable parameters and robust simulations of both streamflow and soil moisture states.

**Plain Language Summary** Soil moisture plays a critical role in understanding water movement in the environment, influencing processes like plant growth and groundwater recharge. This study compares two ways of modeling water movement in the soil at a regional scale. The first scheme uses a complex model called the Richards equation (RE), which simulates both upward and downward water movement through the soil. The second scheme, called infiltration capacity (IC), is simpler and only allows water to move downward. We found that the RE scheme is better at predicting soil moisture content and performs similarly to the IC-based model in simulating river discharge, though it requires more computational time. The IC scheme is faster and can still accurately predict soil moisture anomalies, making it suitable for applications that need rapid results, such as flood forecasting and drought monitoring.

## 1. Introduction

Soil moisture (SM) infiltration is a critical process in hydrologic modeling, as it directly affects streamflow generation, groundwater recharge, evapotranspiration, and even biogeochemical processes such as nutrient transport (Koster et al., 2004; Nearing et al., 2016; Seneviratne et al., 2010; Vereecken et al., 2019). Hydrologists working with mesoscale models, such as land surface models (LSMs) and global hydrologic models (GHMs) which operate at grid scales ranging from 1 to 10,000 km<fnref n="2" />, implement SM infiltration using contrasting approaches. One class of models applies the one-dimensional Richards equation (1D RE), which simulates two-way vertical fluxes driven by both pressure head gradients and gravitational drainage, as used in models like MAT-SIRO (Koirala et al., 2014), LEAF2-Hydro (Miguez-Macho et al., 2007), and VIC-ground (Liang & Xie, 2003). In contrast, other models adopt the infiltration capacity (IC) approach, which simplifies infiltration as a one-way downward flux, neglecting capillary rise or upward movement, as seen in models like PCR-GLOBWB (Beek & Bierkens, 2009), HBV (Bergström, 1995), and VIC-2L (Liang et al., 1994).

<footnote n="1">¹Department Environmental Informatics, UFZ - Helmholtz Centre for Environmental Research GmbH, Leipzig, Germany,</footnote>

<footnote n="2">²University of Potsdam, Institute of Environmental Science and Geography, Potsdam, Germany, ³Department</footnote>

<footnote>Computational Hydrosystems, UFZ - Helmholtz Centre for Environmental Research GmbH, Leipzig, Germany, ⁴Faculty of Environmental Sciences, Czech University of Life Sciences, Prague, Czech Republic, ⁵Université de Lorraine, AgroParisTech, INRAE, Nancy, France</footnote>
**Methodology:** A. Kholis, O. Rakovec, L. Samaniego

**Project administration:** T. Kalbacher, L. Samaniego

**Resources:** T. Kalbacher, L. Samaniego

**Software:** A. Kholis, M. Cuntz, L. Samaniego

**Supervision:** T. Kalbacher, L. Samaniego

**Validation:** A. Kholis, O. Rakovec, F. Boeing

**Visualization:** A. Kholis

**Writing – original draft:** A. Kholis

**Writing – review & editing:** A. Kholis, T. Kalbacher, O. Rakovec, F. Boeing, M. Cuntz, L. Samaniego

One of the key advantages of using the 1D RE scheme is its capability to capture vertical interactions between groundwater, soil moisture, and the atmosphere. These exchanges are critical when feedback between different Earth system components are a primary concern, particularly in LSMs and GHMs that continue to advance their hydrological process representations (Clark et al., 2015; Döll et al., 2016; Gleeson et al., 2021). However, implementing the 1D RE at the mesoscale poses challenges due to scale mismatch, specifically, how to relate the equation valid at point scale into appropriate parameters and process representations at coarser spatial resolutions (Blöschl et al., 2019; Wood et al., 2011). The RE was originally formulated and derived under the assumption of homogeneity and local equilibrium, which typically hold only at fine spatial scales, for example, less than 10 m (Or et al., 2015). Additionally, this equation possesses a significant degree of nonlinearity resulting from the pressure dependency of the parameters of soil moisture and hydraulic conductivity (Farthing & Ogden, 2017; Paniconi et al., 1991).

These limitations—stemming from the local equilibrium assumption and the inherent nonlinearity—often lead to biased results, as the RE implementation does not explicitly account for subgrid variability at mesoscale resolutions (Lehmann et al., 2008; Lehmann & Or, 2009; Vogel & Ippisch, 2008). This challenge in developing appropriate scale-dependent subgrid parameterizations for RE-based models may partly explain their underperformance compared to calibrated conceptual models (e.g., Beck et al., 2016; Haughton et al., 2016; Hou et al., 2023; Krysanova et al., 2020), although this is also influenced by uncertainties in model structures, process representations, and parameterizations beyond infiltration alone (Hou et al., 2023). A widely proposed solution to address this challenge is to identify “effective parameters,” defined as parameters that represent the heterogeneity of the grid, through a parameter estimation process (Beven, 1989; Kirchner, 2006).

However, the application of 1D RE in mesoscale models rarely involves such calibration process, primarily due to: (a) the prevailing confidence in the general applicability of physics-based equations, which leads to the use of fixed, hard-coded parameter values or constant pedo-transfer function (PTF) coefficients without calibration (Cuntz et al., 2016; Göhler et al., 2013; Mendoza et al., 2015), and (b) the lack of widely accepted regionalization techniques to address subgrid variability and nonlinearity, as well as to reduce the dimensionality of parameter fields to avoid overparameterization (Beck et al., 2016; Samaniego et al., 2017). In contrast, the IC scheme, commonly applied in conceptual models, is parsimonious and computationally efficient, making it more practical for parameter calibration and regionalization (Klemeš, 1983; Mendoza et al., 2015). Therefore, calibrating the 1D RE to identify its effective parameters represents a promising step toward improving its performance in mesoscale modeling. Developing improved parameterization strategies for the RE is essential for enhancing LSMs and, in turn, for better understanding the impacts of climate change on the Earth system (Samaniego, 2025).

One potential method to obtain effective parameters for the 1D RE is the Multiscale Parameter Regionalization (MPR) technique, as proposed by Samaniego et al. (2010). MPR follows the principle of “estimate the parameters first, then average, but not necessarily arithmetically,” where the parameterization is based on established and well-tested transfer functions (TFs) (e.g., PTFs) to estimate soil hydraulic properties. This approach addresses subgrid variability, nonlinearity, and overparameterization simultaneously, ensuring consistent fluxes and state variables when the parameters are transferred across different spatial scales and basins (Kumar, Samaniego, & Attinger, 2013; Samaniego et al., 2017). The regional transferability of MPR has been demonstrated by Kumar, Livneh, and Samaniego (2013), who successfully applied parameters derived from a single Central European basin to 80 candidate basins in the United States. Moreover, MPR has been applied for multivariate parameter estimation and evaluation, combining various observation types such as streamflow, soil moisture, and evaporation (Rakovec et al., 2019; Rakovec, Kumar, Attinger, & Samaniego, 2016; Zink et al., 2018).

An early attempt to regionalize 1D RE parameters using the MPR was carried out in the VIC model (Mizukami et al., 2017; Newman et al., 2017), which showed improvements in the spatial distribution of model parameters and modest improvements in simulated streamflow. However, that study focused only on saturated hydraulic conductivity, leaving other important RE parameters, such as saturated matric potential and pore size distribution index, unchanged. Additionally, the same PTFs were retained without further evaluation. In fact, the selection of PTFs for the RE parameters represents a pivotal stage, as the choice of PTFs has been shown to significantly influence soil moisture estimates and, consequently, the simulated water fluxes (Soylu et al., 2011; Weihermüller et al., 2021). Furthermore, to the best of the authors' knowledge, no previous studies have focused on producing transferable RE parameters that work consistently across different locations and spatial scales.
This study has two main objectives: (a) to demonstrate that the 1D RE-based formulation can be made scalable and transferable using MPR technique, and (b) to compare the RE-based model (mHM-RE) against the IC-based scheme (mHM-IC) to assess whether such physically based approach can perform competitively at the mesoscale. This research serves as an important step in laying the foundation for scalable, physics-based hydrologic modeling—particularly toward coupling with groundwater models for continental-to global-scale applications within the mHM framework.

To address these objectives, the 1D RE parameterization was carried out by testing three different PTFs. Parameter transferability was evaluated through cross-validation across multiple basins and spatial scales. Model performance was assessed using streamflow observations from 201 basins and soil moisture measurements from 46 sites (0–25 cm depth) and 42 sites (25–60 cm and 0–60 cm depths) across Germany. All simulations were performed using the mHM, which natively incorporates MPR and typically represents soil infiltration using the IC approach (Samaniego et al., 2010). The 1D RE was implemented in mHM following the Soil-Litter-Iso (SLI) model developed by Haverd and Cuntz (2010).

## 2. Modeling Schemes and Parameterization

### 2.1. The mHM-IC and mHM-RE Models

The mesoscale Hydrologic Model (mHM; Samaniego et al., 2010; Kumar, Samaniego, & Attinger, 2013) was originally a bucket-type, process-based distributed model designed for applications at the mesoscale basin level. This model simulates key hydrologic processes, including ET, canopy interception, the accumulation and melting of snow, infiltration, soil moisture dynamics, runoff generation, and routing. For detailed overviews and parameterizations of the model, refer to Kumar (2010) and Samaniego et al. (2010). The advantage of using mHM compared to other models is that it already provides a platform for testing various settings for all state variables and key flows. The implementation of the Richards-based module in mHM was straightforward, as all other components of the model were kept constant. This consistency is critical for enabling a statistically sound comparison. Such controlled conditions are rarely achievable in previous intercomparison projects for example, ISI-MIP (Warszawski et al., 2014), DMIP (Smith et al., 2004, 2012), and PILPS (Henderson-Sellers et al., 1993, 1995), which often involve multiple preprocessing steps such as data aggregation, hard-coded TF parameters—commonly found in LSMs like NOAH-MP and CLM (Mendoza et al., 2015)—and varying numerical integration schemes. In contrast, the present experimental setup avoids these limitations, allowing for the first time an isolated assessment of the two soil process modules and their respective parameterizations.

The mHM is driven by a set of input variables including precipitation (*P*), temperature (*T*), and potential evapotranspiration (*PET*). The processes associated with this study happen individually on the single grid cells and are illustrated in Figure 1. Total evapotranspiration (*ET*) is composed of three components: canopy evaporation (*E*<sub>C</sub>), sealed area evaporation (*E*<sub>D</sub>), and soil layer evapotranspiration (*E*<sub>i</sub>, *i* = 1–6 in this study). Effective precipitation (*P*<sub>eff</sub>), encompassing throughfall as well as snow melting, is distributed to either soil layers or sealed areas depending on the land cover. In case of impervious surfaces, *P*<sub>eff</sub> is directed toward the sealed area bucket and subsequently treated as direct runoff (*Q*<sub>D</sub>) and *E*<sub>D</sub>. For pervious surfaces, *P*<sub>eff</sub> infiltrates into the first soil layer, and following that, SM dynamics across multiple layers (1 to *n*) are calculated. The vertical flux originating from the last soil layer (*q*<sub>n</sub>) is directed toward the subsurface reservoir, where fast and slow interflow (*Q*<sub>If</sub> and *Q*<sub>Is</sub>) as well as groundwater recharge are generated. The groundwater reservoir is fed by recharge, and the calculation of baseflow (*Q*<sub>B</sub>) from this compartment depends on the storage level. The runoff components *Q*<sub>D</sub>, *Q*<sub>Is</sub>, *Q*<sub>If</sub>, and *Q*<sub>B</sub> within a given cell are first summed to obtain the total runoff (*Q*), which is then routed toward the basin outlet (Thober et al., 2019).

In the original version of mHM (Samaniego et al., 2010), soil moisture (SM) transmission is simulated using the IC approach, where the vertical flux (*q* − IC in Figure 1) permits water movement only in the downward direction. These vertical fluxes are activated only during effective precipitation events. In this study, the one-dimensional RE (1D RE) is introduced into mHM to simulate the bidirectional dynamics of soil moisture (SM), represented as *q*-RE in Figure 1. This is achieved by integrating a module from the Soil-Litter-Iso (SLI) model, developed by Haverd and Cuntz (2010) and Cuntz and Haverd (2018). SLI is a one-dimensional model that simulates the coupled transport of heat, water, and stable isotopes through soil and litter layers. Its capabilities allow to trace water flow through soil and snow, partition terrestrial water, and represent permafrost processes
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Diagram of the main processes of mHM related to soil moisture dynamics.</figcaption></figure>

such as freezing and thawing are well established (Cuntz & Haverd, 2018). SLI has also been incorporated into the LSM CABLE (Haverd et al., 2013, 2016). In this study, however, only the components of SLI related to SM dynamics are retained, while processes involving heat and isotopes are excluded. For clarity, the mHM versions using the IC and RE approaches are hereafter referred to as mHM-IC and mHM-RE, respectively. Parameterizations for both model configurations are provided in Section 2.3.

## 2.2. Multiscale Parameter Regionalization (MPR)

The MPR procedure involves two main steps. First, model-effective parameters—such as saturated hydraulic conductivity and residual water content—are derived at the native data scale ($l_0$) using predefined TFs. Second, these parameters are upscaled to the target model resolution ($l_1$). Both the TFs and the associated upscaling operators must be explicitly defined prior to MPR implementation. Each TF consists of three core elements: (a) effective parameters that are subject to upscaling, (b) coefficients—termed global parameters—that are calibrated and transferable across basins, and (c) physiographic predictors such as soil texture and slope. The current MPR
framework supports five types of upscaling operators: arithmetic mean, harmonic mean, geometric mean, maximum difference, and majority rule. A comprehensive description of the MPR methodology is provided in Samaniego et al. (2010) and Schweppe et al. (2022).

In MPR, calibration is applied on the coefficients of the TFs (γ), also referred to as global parameters, which are associated with spatially distributed basin predictors—rather than with the grid-specific effective parameters. This approach substantially reduces the degrees of freedom in parameter optimization compared to non-regionalized calibration methods (Kumar, Samaniego, & Attinger, 2013; Samaniego et al., 2010). By calibrating the TF coefficients, effective model parameters can be estimated at the native data resolution ($l_0$), and subsequently upscaled to any desired target resolution ($l_1$).

MPR requires predefined TFs and upscaling operators in order to estimate effective parameters. To find a robust TF for the one-dimensional Richards equation (1D RE), this study evaluates three different pedo-transfer functions (PTFs). These PTFs are used to compute five newly introduced effective parameters within the mHM-RE model: saturated hydraulic conductivity ($K_s$), saturated water content or porosity ($\theta_s$), residual water content ($\theta_r$), air entry potential ($h_s$), and pore size distribution index ($\lambda$). The corresponding TFs are summarized in Table 1 and described in detail in Section 2.3. The choice of upscaling operators was guided by analogies from soil physics, where the physical units of the parameters often inform the most appropriate aggregation method. Following the recommendations of Kumar, Samaniego, and Attinger (2013), the harmonic mean was applied to all newly implemented parameters in this study.

## 2.3. Transfer Functions of Soil Moisture Dynamics

The continuum equation for SM dynamics in an unsaturated layer can be expressed as:

$$ \frac{\partial \theta}{\partial t} = \frac{\partial q}{\partial z} - E \qquad (1) $$

where $\theta$ is volumetric SM (m$^3$ m$^{-3}$), $t$ is time (s), $z$ is soil depth positive downward (m), $q$ is vertical flux (m s$^{-1}$), and $E$ is soil evapotranspiration (s$^{-1}$).

In order to compute $q$, the IC approach applies the concept of potential infiltration rates (e.g., Horton, 1941) for not only the surface soil but for all soil layers. This formulation neglects the diffusivity of vertical flux, thereby restricting the transport of SM to downward movement exclusively. $q$ in IC can originate from the flux of the layer above as in mHM (Samaniego et al., 2010) and HBV (Bergström, 1995), from saturated hydraulic conductivity ($K_s$) as in VIC (Liang et al., 1994), or from unsaturated hydraulic conductivity ($K$) as in PCR-GLOBWB (Beek & Bierkens, 2009). The variable $q$ in mHM-IC is represented mathematically as follows:

$$ q_i = q_{i-1} \left[ 1 - \left( \frac{\theta}{\theta_s} \right)^{\lambda} \right] \qquad (2) $$

where $\theta_s$ is the saturated SM content (m$^3$ m$^{-3}$), $\lambda$ denotes the pore size distribution index (-), and $i$ represents the soil layers from 1 to $n$. The vertical flux feeding the first soil layer ($q_0$) is $P_{eff}$ as discussed before. Equation 2 contains two unknown parameters: $\theta_s$ and $\lambda$. The value of $\theta_s$ is calculated using a PTF from Zacharias and Wessolek (2007) and $\lambda$ is parameterized as a function of the soil bulk density ($D_b$) in mHM-IC (see Table 1). Equation 1 is then solved by following discretization, applied to each soil layer:

$$ \theta_i^k = \theta_i^{k-1} + q_i^k - E_i^k \qquad (3) $$

where $k$ is time step index.

In the 1D RE, the vertical flux of liquid water $q$ is the sum of transport due to a gradient in pressure head $h$ and transport by gravitational drainage. The SLI model (Haverd & Cuntz, 2010), and therefore mHM-RE, applies the Kirchhoff transform to linearize 1D RE using the matric flux potential (MFP, $\Phi$):
<table><caption>**Table 1**

Soil Moisture Relations and Transfer Functions for mHM-IC and the 1D Richards Equation Implementations in mHM-RE</caption><thead><tr><th>Model</th><th>Moisture relations</th><th>Transfer functions</th></tr></thead><tbody><tr><td>mHM-IC</td><td></td><td>Zacharias and Wessolek (2007)<br>θ<sub>s</sub> = {θ<sub>s</sub> = γ<sub>1</sub> + γ<sub>2</sub>C<sub>l</sub> + γ<sub>3</sub>D<sub>b</sub> if S<sub>a</sub> &lt; 66.7<br>θ<sub>s</sub> = γ<sub>4</sub> + γ<sub>5</sub>C<sub>l</sub> + γ<sub>6</sub>D<sub>b</sub> otherwise<br>λ = γ<sub>7</sub>D<sub>b</sub></td></tr><tr><td>mHM-RE 1</td><td>Clapp and Hornberger (1978)<br>S = θ/θ<sub>s</sub> = (h/h<sub>s</sub>)<sup>-λ</sup><br>K = K<sub>s</sub>S<sup>2+3</sup></td><td>Cosby et al. (1984)<br>h<sub>s</sub> = 10<sup>(γ<sub>1</sub> - γ<sub>2</sub>S<sub>a</sub> + γ<sub>3</sub>S<sub>s</sub>)</sup><br>K<sub>s</sub> = 10<sup>(γ<sub>4</sub> + γ<sub>5</sub>S<sub>a</sub> - γ<sub>6</sub>C<sub>l</sub>)</sup><br>λ = γ<sub>7</sub> - γ<sub>8</sub>S<sub>a</sub> + γ<sub>9</sub>C<sub>l</sub><br>θ<sub>s</sub> = γ<sub>10</sub> - γ<sub>11</sub>S<sub>a</sub> - γ<sub>12</sub>C<sub>l</sub></td></tr><tr><td>mHM-RE 2</td><td>Clapp and Hornberger (1978)<br>S = θ/θ<sub>s</sub> = (h/h<sub>s</sub>)<sup>-λ</sup><br>K = K<sub>s</sub>S<sup>η</sup></td><td>Campbell (1985)<br>h<sub>s</sub> = h<sub>e</sub>(D<sub>b</sub>/1.3)<br>K<sub>s</sub> = γ<sub>1</sub> exp(γ<sub>2</sub>S<sub>s</sub> + γ<sub>3</sub>C<sub>l</sub>)<br>λ = σ<sub>g</sub> - h<sub>e</sub><br>θ<sub>s</sub> = γ<sub>4</sub>(1 - D<sub>b</sub>/γ<sub>5</sub>)</td></tr><tr><td>mHM-RE 3</td><td>Brooks and Corey (1964)<br>S = θ - θ<sub>r</sub>/θ<sub>s</sub> - θ<sub>r</sub> = (h/h<sub>s</sub>)<sup>-λ</sup><br>K = K<sub>s</sub>S<sup>η</sup></td><td>where:<br>h<sub>e</sub> = -1/√d<sub>g</sub><br>σ<sub>g</sub> = exp[-γ<sub>6</sub>S<sub>s</sub> + γ<sub>7</sub>C<sub>l</sub> - ln<sup>2</sup>(d<sub>g</sub>)]<br>d<sub>g</sub> = exp[-γ<sub>8</sub> - γ<sub>9</sub>S<sub>s</sub> - γ<sub>10</sub>C<sub>l</sub>]</td></tr><tr><td></td><td></td><td>Rawls and Brakensiek (1989)<br>h<sub>s</sub> = exp[γ<sub>1</sub> + γ<sub>2</sub>C<sub>l</sub> - γ<sub>3</sub>θ<sub>s</sub> - γ<sub>4</sub>C<sub>l</sub><sup>2</sup> - γ<sub>5</sub>S<sub>a</sub>θ<sub>s</sub> - γ<sub>6</sub>C<sub>l</sub>θ<sub>s</sub> +<br>γ<sub>7</sub>S<sub>a</sub><sup>2</sup>θ<sub>s</sub> - γ<sub>8</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> - γ<sub>9</sub>S<sub>a</sub><sup>2</sup>C<sub>l</sub> + γ<sub>10</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub> - γ<sub>11</sub>S<sub>a</sub><sup>2</sup>θ<sub>s</sub> +<br>γ<sub>12</sub>C<sub>l</sub><sup>2</sup>S<sub>a</sub> + γ<sub>13</sub>θ<sub>s</sub><sup>2</sup>C<sub>l</sub>]<br>K<sub>s</sub> = exp[γ<sub>14</sub>θ<sub>s</sub> - γ<sub>15</sub> - γ<sub>16</sub>C<sub>l</sub> + γ<sub>17</sub>S<sub>a</sub><sup>2</sup> - γ<sub>18</sub>C<sub>l</sub><sup>2</sup> - γ<sub>19</sub>θ<sub>s</sub><sup>2</sup> +<br>γ<sub>20</sub>S<sub>a</sub>θ<sub>s</sub> - γ<sub>21</sub>S<sub>a</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> - γ<sub>22</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> + γ<sub>23</sub>S<sub>a</sub><sup>2</sup>C<sub>l</sub> + γ<sub>24</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub> +<br>γ<sub>25</sub>S<sub>a</sub><sup>2</sup>θ<sub>s</sub> - γ<sub>26</sub>C<sub>l</sub><sup>2</sup>S<sub>a</sub>]<br>λ = exp[-γ<sub>27</sub> + γ<sub>28</sub>S<sub>a</sub> - γ<sub>29</sub>θ<sub>s</sub> - γ<sub>30</sub>S<sub>a</sub><sup>2</sup> - γ<sub>31</sub>C<sub>l</sub><sup>2</sup> + γ<sub>32</sub>θ<sub>s</sub><sup>2</sup> -<br>γ<sub>33</sub>S<sub>a</sub>θ<sub>s</sub> + γ<sub>34</sub>S<sub>a</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> - γ<sub>35</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> - γ<sub>36</sub>S<sub>a</sub><sup>2</sup>C<sub>l</sub><sup>2</sup> +<br>γ<sub>37</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub> - γ<sub>38</sub>C<sub>l</sub>θ<sub>s</sub><sup>2</sup>]<br>θ<sub>s</sub> = γ<sub>39</sub>(1 - D<sub>b</sub>/γ<sub>40</sub>)<br>θ<sub>r</sub> = -γ<sub>41</sub> + γ<sub>42</sub>S<sub>a</sub> + γ<sub>43</sub>C<sub>l</sub> + γ<sub>44</sub>θ<sub>s</sub> - γ<sub>45</sub>C<sub>l</sub><sup>2</sup> - γ<sub>46</sub>S<sub>a</sub>θ<sub>s</sub> - γ<sub>47</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub><sup>2</sup> +<br>γ<sub>48</sub>C<sub>l</sub><sup>2</sup>θ<sub>s</sub> - γ<sub>49</sub>C<sub>l</sub>θ<sub>s</sub><sup>2</sup></td></tr></tbody></table>

Note. γ is global (calibrated) parameter (-), h_s is air entry potential (m), K_s is saturated hydraulic conductivity (m s⁻¹), λ is the pore distribution index, θ_s is saturated soil moisture (m³ m⁻³), θ_r is residual soil moisture (m³ m⁻³), S is degree of effective saturation (-), S_a is percentage of sand, S_l is percentage of silt, C_l is percentage of clay, D_b is bulk density (g cm⁻³), h_e (m) is air entry value at the standard D_b = 1.3 (g cm⁻³), d_g is the geometric mean particle diameter (-), and σ_g is the standard deviation of d_g (-).
$$q = -K(h)\left(\frac{\partial h}{\partial z} - 1\right) = -\frac{\partial \Phi}{\partial z} + K(h). \qquad (4)$$

The MFP (Φ), as proposed by Gardner (1958) and Campbell (1985) for the linearization of 1D RE, is expressed as:

$$\Phi = \int_{-\infty}^{h} K(\bar{h}) d\bar{h} = \begin{cases} \frac{K_s h_s}{1 - \lambda\eta} \left(\frac{h}{h_s}\right)^{1-\lambda\eta} & \text{if } h < h_s \\ \frac{K_s h_s}{1 - \lambda\eta} \left(\frac{h}{h_s}\right) & \text{if } h \ge h_s \end{cases} \qquad (5)$$

where η is related to the pore-size distribution parameter λ by η = 2/λ + 3 (Campbell, 1974). Using Φ instead of h in Equation 4 improves the numerical stability and accuracy of flux q calculations, especially when applied to models with coarse vertical discretization, such as LSMs and GHMs. This is particularly important since these models often represent soil layers at much larger thicknesses than typically used in fine-scale soil physics studies. Furthermore, the mHM-RE uses the soil water retention curve and hydraulic conductivity relation of Brooks and Corey (1964), as well as the simplification of Clapp and Hornberger (1978), which assumes a residual water content (θr) of zero:

$$S = \frac{\theta - \theta_r}{\theta_s - \theta_r} = \left(\frac{h}{h_s}\right)^{-\lambda} \qquad (6)$$

$$K = K_s S^{\eta} \qquad (7)$$

where S is degree of effective saturation (-).

The SLI model uses Ross's fast numerical solution (Ross, 2003) to avoid the large number of iterations required to achieve convergence, as is the case with standard algorithms (e.g., HYDRUS-1D model) for solving 1D RE, which assume that the flux is driven by the gradient of the hydrologic head. Ross's solution, on the contrary uses the MFP as driving force for flow and linearizing q in terms of ΔS. The general discretization to solve Equation 1 for the mHM-RE model is given as:

$$\frac{\Delta z_i (\theta_{r,i} + (\theta_{s,i} - \theta_{r,i}) \Delta S_i)}{\Delta t} = q_{i-1}^{\sigma} - q_i^{\sigma} - E_i \qquad (8)$$

where fluxes $q_{i=1,...,n-1}$ are to be evaluated at fraction σ through the time step. Term $q_i^σ$ is defined as:

$$q_i^\sigma = q_i + \sigma \left( \frac{\partial q_i}{\partial S_i} \Delta S_i + \frac{\partial q_i}{\partial S_{i+1}} \Delta S_{i+1} \right). \qquad (9)$$

Substituting Equation 9 into Equation 8 leads to a tridiagonal set of n linear equations:

$$\begin{align} a_i \Delta S_{i-1} + b_i \Delta S_i + c_i \Delta S_{i+1} &= d_i \nonumber \\ i &= 1, \dots, n \nonumber \\ a_i &= \frac{\partial q_{i-1}}{\partial S_{i-1}} \nonumber \\ b_i &= \frac{\partial q_{i-1}}{\partial S_i} - \frac{\Delta z_i (\theta_{s,i} - \theta_{r,i})}{\sigma \Delta t} \nonumber \\ c_i &= -\frac{\partial q_i}{\partial S_{i+1}}, \quad d_i = -\frac{q_{i-1} - q_i}{\sigma} + \frac{E_i \Delta z_i}{\sigma}. \tag{10} \end{align}$$

This tridiagonal system is efficiently solved using the Thomas algorithm (Bittelli et al., 2015; Campbell, 1985). Equation 10 applies only under unsaturated conditions. When the soil layer becomes saturated, ΔS<sub>i</sub> reduces to
zero, and consequently, $\Delta\Phi_i$ is used instead by replacing $S_i$ with $\Phi_i$ in the computation. To improve numerical stability and accuracy, the weighting factor $\sigma$ is set to 1.0 for saturated layers and 0.5 for unsaturated layers (Haverd & Cuntz, 2010; Ross, 2003). A detailed description of the discretization and numerical solution is provided in Ross (2003) and Haverd and Cuntz (2010).

The selection of PTFs for estimating unknown parameters ($K_s$, $\theta_s$, $\lambda$, $h_s$, and $\theta_r$) in the mHM-RE is guided by two main criteria. First, the PTFs must be formulated to estimate the parameters of the model by Brooks and Corey (1964) or its simplified variant by Clapp and Hornberger (1978), without requiring complex transformations—such as converting the retention curve (Equation 6) into logarithmic form or dividing it into distinct dry and wet regimes (e.g., Guber et al., 2006)—since the Ross solver used in this study accepts only the original formulation of the retention curve. Second, the PTFs must be expressed in continuous functional forms (e.g., regression-based) enabling their use as TFs within MPR, as opposed to class-based PTFs that rely on lookup tables. Based on these criteria, and informed by commonly used PTFs in the literature (e.g., Dai et al., 2019; Guber et al., 2006; Weihermüller et al., 2021; Zhang et al., 2020), this study evaluates three candidate PTFs: those proposed by Cosby et al. (1984), Campbell (1985), and Rawls and Brakensiek (1989), as summarized in Table 1.

These three PTFs differ in both complexity and the number of samples used in their development (Table 1). mHM-RE 1 and mHM-RE 2 omit the residual water content $\theta_r$. The simplest setup, mHM-RE 1 employs the PTFs of Cosby et al. (1984), which relies solely on soil texture as input and introduce 12 global parameters. mHM-RE 2 uses the PTFs of Campbell (1985), which incorporate both soil texture and bulk density ($D_b$), resulting in 10 global parameters. The most complex configuration, mHM-RE 3, is based on the work of Rawls and Brakensiek (1989), includes $\theta_r$, and introduces 49 global parameters (Table 1).

The PTFs also vary in their empirical foundation: Cosby et al. (1984) and Rawls and Brakensiek (1989) developed their functions using 1,488 and 5,371 soil samples, respectively, both from across the United States. In contrast, the PTF by Campbell (1985) is based on only six samples collected near Washington, United States (Campbell & Shiozawa, 1992). This variation in model complexity and data support underscores the importance of evaluating the influence of these PTFs on model performance.

The soil evapotranspiration $E_i$ ($i = 1, \dots, 6$ soil layers used in this study), used in Equation 1 for all model variants, is computed as a function of potential evapotranspiration (PET), canopy evaporation ($E_C$), and soil-stress factor ($F$):

$$E_i = \begin{cases} PET - E_C & \text{if } i = 1 \\ PET - E_C - \sum_{j=1}^{i-1} E_j & \text{if } i > 1 \\ 0 & \text{if } \sum_{j=1}^{i-1} E_j = 0 \end{cases} \quad (11)$$

$$E_i = E_i \cdot F_i \quad (12)$$

where $i$ represents soil layer and $j$ denotes the index of the preceding layers. The $F$ is a soil-stress factor that controls the distribution of the remaining potential evapotranspiration ($PET - E_C$) across each soil layer as function of water content ($\theta$), soil wilting point ($\theta_{wp}$), field capacity ($\theta_{fc}$), and fraction of roots ($f_r$), following the Feddes reduction function (Feddes et al., 1976):

$$F = \begin{cases} f_r & \text{if } \theta \ge \theta_{fc} \\ f_r \frac{\theta - \theta_{wp}}{\theta_{fc} - \theta_{wp}} & \text{if } \theta < \theta_{fc} \\ 0 & \text{if } \theta < \theta_{wp} \end{cases} \quad (13)$$

$$f_r = 1 - \gamma_l^z \quad (14)$$
where z is the soil depth (m, positive downward) and $γ_l$ is a land cover-dependent calibration parameter, with specific values for forest ($γ_1$), impervious cover ($γ_2$), and pervious cover ($γ_3$), as listed in Table S3 ($γ$ soil evapotranspiration) in Supporting Information S1. The cumulative subtraction term in Equation 11 causes the first or shallower soil layer to typically have a higher $E_l$ value than the deeper layers, unless F increases drastically, which is unlikely since $F_l$ is constrained by root fraction $f_r$. Note that all models in this study have similar PET values, as this variable is calculated based on $T_{min}$ and $T_{max}$ input data.

The parameters $θ_{fc}$ and $θ_{wp}$ in Equation 13 are estimated using the PTFs from Twarakavi et al. (2009):

$$θ_{fc} = θ_s \cdot \exp(γ_1 \cdot (γ_2 + \log(K_s)) \cdot n) \quad (15)$$

$$θ_{wp} = \frac{θ_s}{\exp(m \cdot \log(1 + \exp(n \cdot \log(γ_1 \cdot α)))))} \quad (16)$$

where n, m, and α are Genuchten (1980) empirical parameters, computed in mHM following the PTFs by Zacharias and Wessolek (2007):

$$m = 1 - \frac{1}{n} \quad (17)$$

$$n = γ_1 - γ_2 \cdot S_a^{γ_3} + γ_4 \cdot C_l^{γ_5} \quad (18)$$

$$α = \exp(γ_1 + γ_2 \cdot S_a + γ_3 \cdot C_l - γ_4 \cdot D_b) \quad (19)$$

where $γ$ is optimized (global) parameter, $S_a$ is the percentage of sand, $C_l$ is the percentage of clay, and $D_b$ is the soil bulk density (g·cm⁻³).

## 2.4. Study Area and Data Sets

This study was conducted across 201 river basins distributed throughout Germany, with surface areas ranging from $10^1$ to $10^4$ km² (Boeing et al., 2022). These basins span diverse geographical regions, from the northern coastal lowlands to the southern Alpine uplands (Figure 2). The model setup utilized a range of high-resolution data sets, including terrain attributes from the SRTM digital elevation model (Jarvis et al., 2008), land cover from CORINE (EEA, 2020), leaf area index (LAI) from MODIDS (Myneni et al., 2015), soil physical properties from SoilGrids (Poggio et al., 2021), and lithological information from GLIM (Hartmann & Moosdorf, 2012b). Further details are provided in Text S1 and Table S1 in Supporting Information S1.

The spatial resolution of the geomorphologic and vegetation data sets, which represents subgrid variability within the model domain, varies across sources. All such data sets were resampled to a uniform 200 m × 200 m grid, referred to as the native data resolution ($l_0$). The meteorological forcing data—including daily precipitation, mean air temperature ($T_{mean}$), and maximum and minimum daily air temperatures ($T_{max}$ and $T_{min}$)—were obtained from the German Weather Service (DWD) and interpolated to a 1 km × 1 km grid using External Drift Kriging (Samaniego et al., 2011). Potential evapotranspiration (PET) was computed within the model using the Hargreaves method Hargreaves and Samani (1985), which relies on $T_{max}$ and $T_{min}$. Simulations were conducted at multiple target resolutions ($l_1$), specifically 12, 6, 3, and 1 km grids, as detailed in the experimental design section.

Daily observed streamflow data for 201 basins were obtained from the Global Runoff Data Centre (GRDC). In addition, daily absolute soil moisture (SM) observations from 46 sites (0–25 cm depth) and 42 sites (25–60 cm and 0–60 cm depths) across Germany were sourced from TERENO (Zacharias et al., 2024) and FLUXNET (Pastorello et al., 2020), as reported in Boeing et al. (2022). The SM data span varying time periods and were obtained using four methods: cosmic-ray neutron sensing (CRNS), spatially distributed measurements (SDM), single profile measurements (SPM), and lysimeters (LYSI). In this study, FLUXNET sites use only the SPM method, while TERENO sites utilize all four. To ensure data quality, snow cover periods were excluded from the TERENO CRNS observations (Bogena et al., 2022). All available SM measurements were vertically aggregated using weighted averages to match the model's soil layer configuration (described in Section 2.5): 0–25 cm (layers 1–3, L1–L3), 25–60 cm (layer 4, L4), and 0–60 cm (layers 1–4, L1–L4). Additional details on soil moisture data sets are provided in Text S1 and Table S2 in Supporting Information S1.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Map of the study area showing 201 basins and the location of soil moisture measurements from FLUXNET observation network (Pastorello et al., 2020) and four TERENO sites (Zacharias et al., 2024): Central Germany (CG), Rur/Eifel (Rur), Northeast (NE), and the Pre-Alpine Observatory (PAO).</figcaption></figure>

**2.5. Experimental Design and Model Setup**

The mHM-IC and mHM-REs model variants were configured using an identical setup to ensure consistency in the evaluation. This included the same basin physical properties, meteorological forcing, and soil layering scheme, consist of six layers: 0–5 cm (layer 1, L1), 5–15 cm (layer 2, L2), 15–25 cm (layer 3, L3), 25–60 cm (layer 4, L4), 60–100 cm (layer 5, L5), and 100–200 cm (layer 6, L6). Both configurations assumed free drainage as the soil layer's lower boundary condition of the soil column. While the mHM-RE variant supports a range boundary conditions types, including groundwater interactions, these features were not considered in the present study and will be addressed in a follow-up study.

In this study, the effective parameters were computed and optimized at 200 m resolution based on basin predictors and TFs. These parameters were then upscaled to the desired simulation resolutions, following the workflow outlines below:

1. model calibration against streamflow observations was performed at a spatial resolution of 12 km for all 201 basins;

2. cross-basin validation was conducted at the same 12 km resolution
3. the calibrated models at 12 km resolution were subsequently refined to 6, 3, and 1 km resolutions for cross-scale validation across all basins; and

4. point-based soil moisture (SM) validation was carried out at each observation site using a single 200 m grid cell, consistent with the native resolution $l_0$ of the predictor data sets used in this study.

The experimental framework used to assess the performance of mHM-IC and mHM-REs model variants is outlined below:

1. Individual basin calibration was performed against daily observed streamflow over a 30-year period (1989–2018), following a 5-year spin-up, at a spatial resolution of 12 km × 12 km. All global parameters associated with mHM processes (Figure 1) were included in the optimizations, with the only differences between model variants being the parameterization of soil infiltration processes (see Section 2.3). The number of global parameters varied by model configuration: mHM-IC required 49 global parameters, while mHM-RE 1 to 3 involved 51, 48, and 87 parameters, respectively (see Text S2 and Tables S3–S7 in Supporting Information S1). The objective of this single-basin calibration was to allow each model variant to identify its optimal parameter set independently across 201 domains. Calibration was performed using the Dynamically Dimensioned Search algorithm (Tolson & Shoemaker, 2007) with 1 000 iterations. Model performance was evaluated using the Kling-Gupta efficiency (KGE) metric (Gupta et al., 2009) as the objective function. To ensure physical plausibility, the following parameter constraints were imposed during calibration: $\theta_s$ (0.1–0.55 m$^3$ m$^{-3}$), $\theta_r$ (<0.2, m$^3$ m$^{-3}$), $K_s$ (>10$^{-4}$ cm h$^{-1}$), $h_s$ (<200 cm), and $\lambda$ (0.01–1.0). Violations of these bounds were penalized within the objective function. In addition to overall KGE values, the three individual components of KGE were analyzed: Pearson correlation coefficient (r), the bias ratio ($\beta$; mean simulated/mean observed), and the variability ratio ($\alpha$; standard deviation of simulations/observations).

2. To assess the transferability of the model parameters for streamflow prediction, cross-validation was conducted across both basins and spatial scales using the same 30-year period as the calibration phase. For the cross-basin validation, the parameter set calibrated for each basin at 12 km resolution was applied to the remaining 200 basins, also at 12 km resolution. This process also identified the best-performing overall parameter set for each model variant, determined by selecting the one with the highest median KGE value based on the cumulative distribution function of performance metrics across all basins. The best parameter set for each model was subsequently used for soil moisture validation and Germany-wide seamless simulations (see items 4 and 5 below). In addition, cross-scale validation was conducted by applying each basin's optimal 12 km parameter set for simulations at finer resolutions of 6, 3, and 1 km.

3. The optimized RE parameters obtained from step 1 were compared to their original empirical values to assess the plausibility of the calibrated results. Specifically, the post-calibration parameter values were evaluated against those derived from the respective pedo-transfer functions (PTFs)—Cosby et al. (1984) for mHM-RE 1, Campbell (1985) for mHM-RE 2, and Rawls and Brakensiek (1989) for mHM-RE 3—to determine whether the optimization produced physically realistic values.

4. To investigate differences between the model variants, Germany-wide simulations were conducted using seamless parameter fields derived from the cross-basin validation. Variations in fluxes (e.g., runoff, ET, and their individual components) and key parameters (e.g., $\theta_s$ and $K_s$) were analyzed using graphical representations.

5. Daily soil moisture (SM) validation was performed at 46 sites for the 0–25 cm depth and 42 sites for the 25–60 cm and 0–60 cm depths across Germany (Figure 2). The validation covered varying time periods and three soil depth categories: shallow (0–25 cm), deep (25–60 cm), and integrated (0–60 cm) layers (see Text S1 and Table S2 in Supporting Information S1). Model variants were run using a single 200 m × 200 m grid cell ($l_0$ scale) at each observation site to best represent point-scale measurements. The evaluation was based on both volumetric (absolute) SM and SM anomalies. Anomalies were computed by subtracting the moving average (30-day window) from the SM time series to reduce uncertainty associated with short records. Model performance for volumetric SM was assessed using the KGE and its components—Pearson correlation (r), bias ratio ($\beta$), and variability ratio ($\alpha$). Due to large bias and variability in the absolute SM values, the SM anomalies were evaluated using only the correlation coefficient (r).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Individual basin calibrations: cumulative frequency distribution of daily streamflow performance across 201 basins, simulated at a 12 km grid resolution, showing (a) Kling-Gupta efficiency *KGE*, (b) Pearson correlation coefficient *r*, (c) variability *α*, and (d) bias *β*.</figcaption></figure>

# 3. Results and Discussion

## 3.1. Calibration and Cross-Validation

The streamflow performance of the four different models following the individual basin calibration process is presented in Figure 3. The KGE scores of all models exhibit slight variation and demonstrate excellent results, as indicated by median values ranging from 0.85 to 0.90. Only 3 basins (Schwarzbach, Aich, and Wandse) out of 201 had scores below 0.50 in all models. The suboptimal performance in these basins is likely due to human-induced factors, as they are located in densely populated areas, and two dams are operated in the Aich River basin (Dolich et al., 2024). Anthropogenic influence, such as reservoir operations and irrigation practices, were not incorporated into the model setup. For further information on the implementation of reservoir impacts in mesoscale hydrology models, please refer to Shrestha et al. (2024).

In more detail, the mHM-IC configuration shows slightly superior KGE performance compared to the other models variants, as evidenced by its highest median value in Figure 3b, indicating stronger correlation (*r*) with observed streamflow. Following closely are mHM-RE 1 and mHM-RE 2, while mHM-RE 3 shows the weakest *r* score. Meanwhile, all models exhibit comparable variability (*α*) and bias (*β*) scores, as indicated by the dominant vertical line at a value of 1 in the plots (Figures 3c and 3d).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Median Kling-Gupta efficiency values based on streamflow across 201 basins for each model at various grid resolutions. Individual basin calibration conducted at a 12 km (100%) scale, with the optimized parameters transferred to 6 km (50%), 3 km (25%), and 1 km (8.33%) scales.</figcaption></figure>

The global parameters obtained from the streamflow calibration (individual basin calibration) at the 12 km grid scale were subsequently transferred to finer spatial resolutions of 6, 3, and 1 km. Figure 4 shows that all models variants showed a slight decline in performance during cross-scale validation, with the median KGE value decreasing by approximately 0.02—except for the mHM-RE 3 model, which exhibited a slightly larger decline of up to 0.04. This result aligns with findings of Samaniego et al. (2010), where model performance initially improved as the simulation scale was refined to 8 km, but tended to decline modestly at finer resolutions. Despite these minor reductions, the median KGE values remained above 0.8 across models, highlighting the robustness and consistency of both mHM-IC and mHM-RE variants across spatial resolutions.

The optimized parameter sets from each basin were subsequently applied to all other basins at 12 km resolution. This cross-basin validation resulted in 201 cumulative KGE frequency plots per model variant (see Figure S1 in Supporting Information S1). The best-performing runs—identified by the highest median KGE values—are presented in Figure 5, along with three KGE components: Pearson correlation (r), variability ratio (α), and bias ratio (β). Median KGE values are approximately 0.70 for all model variants, except for mHM-RE 2 (green line), which shows a lower value of around 0.60. The relatively weaker performance of mHM-RE 2 is primarily attributed to its reduced correlation *r* compared to the other variants (Figure 5b). Consistent with the results from individual basin calibrations, differences in KGE scores among the model variants are largely driven by the correlation component, whereas bias (α) and variability (β) remain similar across models. This outcome reflects the characteristic of the KGE metric, which—as the calibration objective function—tends to prioritize the reduction of errors in streamflow bias and variability (Mizukami et al., 2019).

In this experiment, cross-basin validation using data from 201 basins resulted in streamflow performance approximately 22% lower than that obtained from individual basin calibration, based on median KGE values. Specifically, the cross-basin validation produced KGE scores ranging between 0.60 and 0.70, whereas individual basin calibration yielded values between 0.85 and 0.90. Previous studies by Kumar, Samaniego, and Attinger (2013) and Samaniego et al. (2010) reported smaller performance reductions of 11% and 15%, respectively. However, those studies were based on significantly fewer basins (45 and 10, respectively), all of which are within the broader set of 201 basins used in this study. Consequently, the earlier experiments likely reflected a narrower range of hydrological and physiographic variability compared to the more diverse set analyzed here. Despite the observed reduction in predictive skill, the median KGE values from cross-basin validation remain above 0.6, indicating that both the mHM-IC and mHM-RE maintain reasonable performance in ungauged basins.

The results further indicate that more complex RE parameterizations tend to produce slightly lower model performance. This is evident in mHM-RE 2 during cross-basin validation, and in mHM-RE 3 during both single-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Cross-basin validation: cumulative frequency distribution of daily streamflow performance across 201 basins, simulated at a 12 km grid resolution, showing (a) Kling-Gupta efficiency *KGE*, (b) Pearson correlation coefficient *r*, (c) variability *α*, and (d) bias *β*.</figcaption></figure>

basin calibration and cross-scale validation. For mHM-RE 2, which relies on both soil texture and bulk density (*D*<sub>b</sub>) to estimate parameters, single-basin calibration based solely on streamflow data may be insufficient to capture the spatial heterogeneity of basin physical characteristics. A multi-basin calibration strategy (e.g., Mizukami et al., 2017; Rakovec et al., 2019) or multivariate parameter estimation using additional observations (e.g., Rakovec, Kumar, Attinger, & Samaniego, 2016; Rakovec, Kumar, Mai, et al., 2016; Zink et al., 2018) may improve the model's transferability across diverse catchments.

In case of the mHM-RE 3, its weaker performance may be attributed to substantially higher number of global parameters (γ) compared to the other model configurations (Table 1). As noted by Samaniego and Bár-dossy (2005), increasing the number of parameters in empirical models tends to elevate model uncertainty. As a result, mHM-RE 3 may require a greater number of calibration iterations than the 1 000 used uniformly across models in this study (see Section 2.5) to achieve optimal performance. Regarding cross-scale performance, the relatively poor results of mHM-RE 3 suggest that alternative upscaling operators may better capture subgrid variability, particularly when using more complex PTFs (Samaniego et al., 2017).

The average runtime performance after calibration, using a single CPU, is summarized in Table 2. In general, runtime increased with basin area, as larger domains require more computational resources. For basin sizes with areas ≤10,000 km², the newly implemented RE modules in mHM exhibited performance comparable to the
<table><caption>**Table 2**
Average Simulation Time (in Minutes) for Forward Runs Using Parameters From Individual Basin Calibrations</caption><thead><tr><th>Area [km²]</th><th>n</th><th>mHM-IC</th><th>mHM-RE 1</th><th>mHM-RE 2</th><th>mHM-RE 3</th></tr></thead><tbody><tr><td>≤500</td><td>99</td><td>6</td><td>7</td><td>4</td><td>7</td></tr><tr><td>500–5,000</td><td>74</td><td>12</td><td>10</td><td>11</td><td>10</td></tr><tr><td>5,000–10,000</td><td>11</td><td>21</td><td>21</td><td>21</td><td>21</td></tr><tr><td>10,000–30,000</td><td>16</td><td>22</td><td>48</td><td>35</td><td>34</td></tr><tr><td>37,809</td><td>1</td><td>26</td><td>119</td><td>165</td><td>58</td></tr></tbody></table>

Note. *n* is the number of simulated basins. All simulations were performed using a single CPU.

default IC scheme. Interestingly, despite solving more complex equations, mHM-RE 2 ran two minutes faster than the mHM-IC for basins areas ≤500 km², and all RE variants slightly outperformed mHM-IC for basins between 500 and 5,000 km². However, substantial performance differences emerged for basin areas >10,000 km². For example, mHM-RE 1 was more than twice as slow as mHM-IC for basins between 10,000 km² and 30,000 km². Moreover, mHM-RE 2 was approximately six times slower than mHM-IC for the largest basin in the data set, which spans 37,809 km².

Previous research by Weihermüller et al. (2021) demonstrated that low saturated hydraulic conductivity $K_s$ values can significantly reduce the numerical performance of the one-dimensional RE (1D RE) using the HYDRUS-1D model. Although the RE solver used in this study differs from that of HYDRUS-1D, similar effects may partly explain the variation in

runtime observed among the mHM-RE variants. In the largest basin (37,809 km² with 263 cells at 12 km × 12 km grid resolution), the average $K_s$ values are 0.18 m/d for mHM-RE 1 and 0.16 m/d for mHM-RE 2—substantially lower than the 3.56 m/d observed in mHM-RE 3. This lower conductivity values contribute to the longer runtimes recorded for mHM-RE 1 (119 min) and mHM-RE 2 (165 min), which are approximately two and three times slower than mHM-RE 3, respectively (see Table 2).

## 3.2. Parameter Comparisons

As outlined in Section 3.1, individual basin calibration can yield a wide range of parameter values that, while varying in physical realism, may still produce models that closely reproduce observed streamflow. In this section, we evaluate the plausibility of the calibrated parameters in the mHM-RE variants by comparing them to the original values derived from three distinct PTFs (Figure 6). This comparison assesses whether the calibration of PTF-based parameters within the mHM framework preserves physically meaningful relationships—that is, whether the resulting parameters remain within the uncertainty ranges reported by the original PTF studies. To generate the plots, the optimized RE parameters ($\lambda$, $h_s$, $\theta_s$, $K_s$, and $\theta_r$) from 39,150 grid cells—covering 201 basins at a 12 km scale—were grouped according to the 11 USDA soil textural classes. For each class, the mean parameter values were calculated and plotted. Note that the figures vary by parameter, unit, and whether a linear or logarithmic scale is applied.

It is important to note that the soil data set used in this study does not include USDA textural classes—namely sand (92% sand), sandy clay (52% sand), silty clay (6% sand), and clay (22% sand)—as shown in Figure 6. The mHM-RE 1 and mHM-RE 3 variants were parameterized using the PTFs from Cosby et al. (1984) and Rawls and Brakensiek (1989), respectively. Both studies provide mean and standard deviation values of the parameters grouped by USDA textural class, enabling direct comparison with the optimized values obtained in this study. In contrast, the PTF by Campbell (1985), used in mHM-RE 2, reports parameter values only for six broader textural categories (Campbell & Shiozawa, 1992). To enable comparability, the original parameter values from Campbell (1985) were reclassified into USDA textural classes using the conversion equations listed in Table 1.

In general, all estimated parameters align well with their original values, as the simulated values fall within either the standard deviation ranges or the 95% prediction intervals. The $\lambda$ parameter which governs the slope of the water retention curve, plays a critical role: higher $\lambda$ values correspond to a steeper decline in water as matric potential decreases, whereas lower values indicate a more gradual reduction. As frequently reported in the literature, $\lambda$ tends to increase with raising sand content—a relationship that holds consistently across all original pedo-transfer functions (PTFs) and mHM-RE model simulations (Figures 6a–6c). It is worth noting that the original data set from Campbell (1985) and Campbell and Shiozawa (1992) showed a negative correlation between $\lambda$ and sand content, likely due to its small sample size of only six soils. However, after reclassifying this data set using the USDA soil texture system, the observed trend reverses and becomes consistent with the other data sets.

The air entry potential ($h_s$) represents a critical suction threshold at which the largest soil pores begin to drain, marking the transition to unsaturated conditions. Finer-textured soils, such as clays, generally exhibit higher $h_s$ values due to their small particle sizes and large specific surface areas, which enhance water retention. All models and original PTFs consistently show a negative correlation between $h_s$ and sand content (Figures 6d–6f).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the mean soil hydraulic properties values for each textural class against sand content between mHM-RE 1, mHM-RE 2, and mHM-RE 3 versus empirical data from Cosby et al. (1984), Campbell (1985), and Rawls and Brakensiek (1989), respectively: (a–c) pore distribution index $\lambda$ [–], (d–f) air entry pressure head $h_s$ [cm], (g–i) porosity $\theta_s$ [m$^3$/m$^3$], (j–l) saturated hydraulic conductivity $K_s$ [cm/h], and (m) residual water content $\theta_r$ [m$^3$/m$^3$]. The blue line and crosses represent results from the model, while the black dots represent the empirical data with the range and 95% prediction intervals. Note that the soil data set used during the calibration of the 201 basins in this study does not include the USDA textural classes of Sand (92% sand), Sandy clay (52% sand), Silty clay (6% sand), and Clay (22% sand).</figcaption></figure>
Similarly, finer soils tend to have higher porosity ($\theta_s$), as their abundance of small particles creates a larger number of small pores. This trend is reflected in all mHM-RE model parameterizations, which show a negative relationship between $\theta_s$ and sand content (Figures 6g–6i). These small pores also impede water movement, leading to lower saturated hydraulic conductivity ($K_s$) in fine-textured soils. Accordingly, both the optimized and original $K_s$ parameterizations display a positive correlation with sand content (Figures 6j–6l). Lastly, the residual water content ($\theta_r$), included only in mHM-RE 3 shows a slight negative correlation with sand content in its optimized $\theta_r$—consistent with the trend in the original PTF (Figure 6m).

The challenges associated with implementing PTFs, particularly regarding extrapolation and scaling, have been extensively discussed by Looy et al. (2017). Extrapolation refers to the application of PTFs beyond the physical domains from which they were derived, where their performance becomes highly dependent on the original data range and the strength of underlying correlations. This limitation is clearly illustrated when different PTFs are applied within the same study area, as in the present research. Here, the use of three distinct PTFs resulted in varying ranges and relationships among optimized parameter values (Figure 6). These discrepancies led to different best-fit parameter sets for the Germany-wide model simulations (Section 3.3) and produced divergent performances in the soil moisture validation (Section 3.4).

A clear example of the extrapolation challenge is provided by the parameterization of the mHM-RE 2 model, which employs the PTF from Campbell (1985). This PTF was originally developed using only six soil samples in the United States. When applied to the 201 basins across Germany in this study, it produced a narrow and low range of saturated hydraulic conductivity ($K_s$ values (0.0005–1.215 cm hr$^{-1}$)). After calibration, this range became even narrower (0.001–0.631 cm hr$^{-1}$), and remained substantially lower than those obtained from the other PTFs and corresponding optimized models.

This behavior is clearly illustrated in the hydraulic conductivity curves $K(\theta)$ (Figures 7h–7n), where the green dashed and solid lines—representing the original Campbell (1985) and mHM-RE 2, respectively—consistently lie below those of the other models across all soil texture classes. The low $K_s$ values in both the original and the calibrated versions are partially compensated by higher optimized values of porosity $\theta_s$, which can be observed in the moisture retention curves $\theta(h)$ under saturated conditions (Figures 7a–7g). For mHM-RE 2, these elevated $\theta_s$ values lead to higher porosity estimates in northern Germany (Figure 10c), resulting in the overestimation of soil moisture, particularly in the deeper soil layer (Figures 11b and 12).

Moreover, the variability of PTF performance across spatial scales has been widely acknowledged (i.e., Pringle et al., 2007). However, Samaniego et al. (2017) demonstrated that by employing an appropriate parameterization framework—specifically, the MPR technique—can yield spatially consistent (i.e., seamless) fluxes and state variables. In this study, we show that the optimized parameters derived from three different PTFs used to implement the RE remain consistent with their original formulations, despite model calibration being performed at a 12 km grid scale, whereas the original PTFs were developed from point-scale data. Cross-scale validation further confirms that applying MPR with different PTFs yields consistent streamflow performance across scales for the mHM model with the newly implemented RE module (Figure 4).

In this study, the calibration of PTF coefficients led to nearly identical moisture retention curves across all mHM-RE variants (dashed lines in Figures 7a–7g), despite notable differences in $K(\theta)$ relationships (dashed lines in Figures 7h–7n) and $\theta_s$ values. For instance, the original Campbell (1985) PTF, characterized by low $K_s$ values, results in a pronounced plateau in $\theta(h)$, particularly for fine-textured soils such as silty clay loam, silt loam, and clay loam (Figures 7e–7g). In contrast, the calibrated mHM-RE 2 variant eliminates this plateau due to adjustments in $\lambda$ and $h_s$ (Figures 6b and 6g), producing steeper $\theta(h)$ slopes for the same texture classes.

This outcome indicates that calibration based on streamflow data tends to harmonize retention properties, promoting consistency in soil moisture dynamics across variants. The relatively invariant shape of the calibrated retention curves likely explains the comparable streamflow performance achieved by all models across the 201 basins (Figure 3). These findings suggest that calibrating PTF coefficients via MPR-based regionalization using may offer a robust alternative to ensemble PTF strategies (Weihermüller et al., 2021; Zhang et al., 2020) for addressing discrepancies arising from extrapolation and scaling effects.

While calibration reduces discrepancies in the shape of the moisture retention curves, the underlying functional structures of the PTFs used in each mHM-RE variant also influence model behavior. Each PTF relies on distinct empirical formulations and combinations of soil properties to estimate the parameters governing the hydraulic
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Moisture retention θ(h) (a–g) and hydraulic conductivity K(θ) (h–n) curves for seven USDA textural classes. The parameters used to generate each curve are derived from Figure 6 and listed in the Tables S8–S10 in Supporting Information S1. Note that the soil data set for the 201 basins analyzed in this study does not include the USDA textural classes of Sand, Sandy clay, Silty clay, and Clay during the calibration processes.</figcaption></figure>

functions—namely, the moisture retention and hydraulic conductivity curves. As a result, even if the calibrated retention curves appear similar, the models may exhibit differing sensitivities to the calibration data. In particular, the observed variations in K(θ) curves and θs values likely stem from the relatively low sensitivity of streamflow-based calibration to parameters such as Ks and θs. These parameters are often weakly constrained when calibration is based solely on discharge, in the absence of complementary observational data such as soil moisture or evapotranspiration.

A key requirement for applying the MPR technique is the availability of predefined TFs to derive effective parameters—such as the three PTFs used in this study to implement the mHM-RE variants. However, uncertainty in selecting appropriate TFs remains a challenge. This issue may be addressed using automatic function discovery methods, such as Function Space Optimization (FSO) proposed by Feigl et al. (2020), which builds on symbolic regression techniques (Klotz et al., 2017). A preliminary implementation of FSO within the MPR framework has been tested in the mHM model for deriving Ks and θfc (Feigl et al., 2022), demonstrating that automatically discovered TFs can achieve performance comparable to standard MPR parameterizations. While the application of FSO to parameterize the mHM's infiltration process presents a promising direction, it lies beyond the scope of this study.

### 3.3. Germany-Wide Simulation

To assess potential model discrepancies, Germany-wide simulations were performed using the best parameter sets identified through the cross-basin validation (Section 3.1). The spatial distribution of total runoff (Q), presented in Figure 8, reveals broadly consistent patterns across all model variants and observations. Higher runoff values concentrated in the southwestern and southeastern regions, while lower values are observed in the east—patterns
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Germany-wide comparison of total runoff (*Q*): (a) mean monthly precipitation, (b) observed mean monthly total runoff *Q*, (c) *Q* from the mHM-IC, and (d–f) *Q* from mHM-RE 1, mHM-RE 2, and mHM-RE 3, respectively. The models were run using the best simulation parameter set of the cross-basin validation at 12 km grid resolution.</figcaption></figure>

that align with findings of Zink et al. (2017). This consistency is expected, given that the KGE scores for the best-performing simulations differ only slightly among models, with median values ranging from 0.60 to 0.70 (Figure 5). The spatial distribution *Q* closely follows the precipitation pattern, as precipitation is the dominant driver of runoff in the study area (Figure 8a).

Although mHM-IC and mHM-RE yield similar patterns and magnitudes of total runoff (*Q*), they differ notably in how this runoff is partitioned within the subsurface (see Figure S2 in Supporting Information S1). In both modeling approaches, *Q* comprises direct runoff (*Q*<sub>D</sub>), fast and slow interflow (*Q*<sub>If</sub> and *Q*<sub>Is</sub>), and baseflow (*Q*<sub>B</sub>). While *Q*<sub>Is</sub> and *Q*<sub>B</sub> represent the dominant components of total runoff, contributions from *Q*<sub>D</sub> and *Q*<sub>If</sub> are comparatively minor. The key distinction lies in the source of dominant runoff: in mHM-IC, *Q* is primarily generated through *Q*<sub>Is</sub> from the subsurface reservoir, whereas in all mHM-RE variants, *Q*<sub>B</sub> originating from the groundwater reservoir is the main contributor. This difference in runoff partitioning is further corroborated by differences in storage dynamics and groundwater recharge, as shown in Figure S3 in Supporting Information S1. Specifically, mHM-IC exhibits substantially higher average subsurface storage than the mHM-RE variants.
Conversely, all mHM-RE models retain more water in the groundwater reservoir due to their higher average groundwater recharge rates. In mHM, groundwater recharge is defined as the flux from subsurface storage to the groundwater reservoir (Figure 1).

The spatial distributions of total evapotranspiration (ET), although not directly validated against observations, exhibit general consistency across the model variants (Figure 9). Higher ET values are concentrated in the southern part of the study area, while lower values occur in the north. In the mHM model, ET includes contributions from sealed surfaces $E_D$, canopy $E_C$, and evaporation/transpiration from soil layers $E_i$ ($i$ = 1–6 in this study). The Germany-wide simulations indicate that a substantial portion of precipitation is intercepted by vegetation, with canopy evaporation ($E_C$) accounting for approximately 30%–40% of total ET (Figures 9i–9l). The remaining ET primarily results from soil evaporation and plant transpiration ($E_i$ in mHM), while sealed areas contribute only minimal—consistent with the small runoff contribution from these areas.

Among variants, mHM-RE 1 displays slightly lower mean monthly ET, mainly due to reduced canopy evaporation ($E_C$) in the northeastern region (Figure 9j). This reduction is attributed to a lower calibrated value of the canopy interception factor ($\gamma$ interception; see Table S3 in Supporting Information S1), which leads to decreased interception and increased throughfall in that area (Figure S4 in Supporting Information S1). These altered surface water dynamics also contribute to slightly higher total runoff (Q) in mHM-RE 1, particularly across northeastern Germany (Figure 8d).

Figures 9m–9p illustrate that mHM-RE variants produce higher soil evapotranspiration $E_i$ compared to mHM-IC. In mHM-RE 1 and mHM-RE 3, these increase is primarily driven by lower canopy evaporation $E_C$, which leaves a larger residual potential evapotranspiration term ($PET - E_C$) available for soil evaporation, as defined in Equation 11. In contrast, mHM-RE 2 exhibits $E_C$ values comparable to those of mHM-IC but still generates greater $E_i$. This can be attributed to the higher values of the soil-stress factor F in mHM-RE 2, as described in Equation 12. The combination of relatively high $E_C$ and elevated $E_i$ results in the highest total evapotranspiration (ET) among the models, particularly in mHM-RE 2. This increased evaporative loss contributes to slightly reduced total runoff (Q), especially in northeastern Germany (Figure 8e).

Figure 10 summarizes the spatial distributions of key parameters governing soil moisture (SM) transmission, including saturated hydraulic conductivity ($K_s$) and porosity ($\theta_s$). All model variants exhibit a similar north–south gradient in $K_s$ (Figures 10a–10d), although the absolute magnitudes vary across models. In contrast, the spatial patterns of $\theta_s$ differ markedly from those of $K_s$ (Figures 10e–10h), with lower $\theta_s$ values generally found in the northern regions. This inverse relationship between $K_s$ and $\theta_s$ is consistent with the soil water retention characteristics described by the soil water characteristic (SWC) equations (Table 1).

The maps further show that mHM-RE 2 exhibits the lowest saturated hydraulic conductivity $K_s$ values, predominantly within the range of 0–0.9 cm/hr. These low values are consistent with the mean optimized $K_s$ values across 201 basins, as shown in Figure 6k, where mHM-RE 2 consistently yields the lowest estimates. These reduced $K_s$ values also contribute to the relatively higher porosity ($\theta_s$) observed in mHM-RE 2, particularly in the northern part of the study area (10g). In comparison, mHM-RE 3 exhibits a broader and generally higher $\theta_s$ range (0.41–0.54) than both mHM-IC (0.28–0.54) and mHM-RE 1 (0.15–0.54), despite having a relatively wider range of $K_s$ values. This behavior is likely attributable to the inclusion of residual water content $\theta_r$ in the SWC formulation of mHM-RE 3. The increase in $\theta_s$ serves to compensate for $\theta_r$, ensuring sufficient available soil moisture within the model.

### 3.4. Soil Moisture (SM) Validation

The performance of the models in simulating absolute soil moisture (SM) was evaluated at 46 sites for the 0–25 cm depth) and at 42 sites (25–60 cm and 0–60 cm depths across Germany (Figure 11). Each model was run using the best parameter sets obtained from the cross-basin validation (see Section 3.1), consistent with the Germany-wide simulations.

Overall, the models exhibit limited skill in reproducing absolute SM values, with median KGE scores below 0.50 and substantial bias ($\beta$) and variability ($\alpha$) across all depth intervals (Figure 11). Nevertheless, the models perform considerably better in capturing temporal dynamics, as indicated by correlation coefficients ($r$) of approximately 0.80 for most cases—except for mHM-IC at the 25–60 cm depth, where performance is noticeably lower (Figures 11d–11f). The models' ability to reproduce SM anomalies is further supported by strong anomaly
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Germany-wide maps of mean monthly total evapotranspiration (ET) and its components in mm/month: (a–d) total evapotranspiration ET, (e–h) evaporation from sealed area $E_D$, (i–l) evaporation from canopy $E_C$, and (m–p) evapotranspiration from soil layer-1 to layer-6 $E_i$ ($i$ = 1–6 in this study).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spatial distribution of model parameters across Germany, averaged across six soil layers (2 m depth): (a–d) saturated hydraulic conductivity $K_s$ and (e–h) saturated volumetric water content $\theta_s$.</figcaption></figure>

correlation values and consistent patterns in the anomaly time series (see Figures S5 and S9–S11 in Supporting Information S1).

Regarding absolute soil moisture (SM) performance, mHM-IC and mHM-RE 1 exhibit the lowest median KGE scores in the shallow soil layer (0–25 cm), while mHM-RE 2 performs worst in the deeper 25–60 cm layer (Figures 11a–11c). Notably, mHM-IC displays a left-skewed KGE distribution across all depth, interval indicating consistently poor performance at several sites. These low KGE values are further reflected in associated performance metrics: mHM-IC shows high variability ($\alpha$; Figures 11g–11i), mHM-RE 1 exhibits low bias values ($\beta$; Figures 11j–11l), and mHM-RE 2 shows elevated bias, particularly in the deeper layers (Figure 11k).

The poor performance of mHM-IC and mHM-RE 2 is primarily linked to overestimation of SM, as evident in the absolute SM time series (Figures S6–S8 in Supporting Information S1). This overestimation is likely driven by elevated porosity $\theta_s$ values in their parameterizations—especially in southern Germany for mHM-IC and northern Germany for mHM-RE 2 (Figure 10). In contrast, mHM-RE 1 underestimate SM, as indicated by its low $\beta$ values and time series patterns (Figure S6 in Supporting Information S1), which can be attributed to comparatively low $\theta_s$ values—dropping to approximately $0.15 \text{ m}^3\cdot\text{m}^{-3}$ in northern Germany (Figure 10). These results highlight a key limitation: calibrating solely against streamflow may be insufficient for accurately estimating $\theta_s$, ultimately compromising the model's ability to simulate absolute SM. This occurs despite the strong performance observed in streamflow simulations, as discussed in Section 3.2.

In addition to its higher variability, mHM-IC exhibits notable limitations in simulating soil moisture (SM) in deeper layers, as reflected by its relatively poor correlation coefficients ($r$) for absolute SM (Figure 11e) and SM anomalies (Figure S5b in Supporting Information S1). These deficiencies can be attributed to structural limitations in the IC method for vertical SM transmission. Specifically, the IC method permits only downward flow of water, and in mHM-IC, soil moisture can be depleted solely via soil evapotranspiration ($E_t$). As a result, deeper soil layers tend to become fully or nearly saturated when $E_t$ is minimal, leading to persistently high moisture levels.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Cumulative frequency distributions based on volumetric soil moisture ($\theta$) over 46 sites (0–25 cm depth) and 42 sites (25–60 cm and 0–60 cm depths) at a daily time step and 200 m grid resolution, obtained from the forward run using the best simulation parameters of the cross-basin validation: (a–c) KGE scores, (d–f) Pearson correlation coefficient *r*, (g–i) variability $\alpha$, and (j–l) bias $\beta$.</figcaption></figure>

It should be noted that the spatial resolution of the mHM variants is on the order of hundreds meters, whereas the observations are obtained at point scale (few centimeters). Potential sources of mismatch also stem from the kilometer-scale forcing data, which may not accurately capture the micro-meteorological conditions affecting soil moisture (SM) observations. Additionally, local soil properties may be poorly represented in the SoilGrids map
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Time series of absolute soil moisture ($\theta$) and daily soil evapotranspiration ($E_i$) at the Aachen site for different soil layers: Layer 1 to 3 (L1–L3: 0–25 cm), Layer 4 (L4: 25–60 cm), Layer 5 (L5: 60–100 cm), and Layer 6 (L6: 100–200 cm). The observations (black) are only present for SM at L1–L3 and L4. For layers L1–L3, $\theta_i$ is depth-averaged while $E_i$ is the sum of $E_1$ through $E_3$.</figcaption></figure>

used in the simulations. Further, there are well-documented uncertainties associated with SM measurements in the literature (e.g., Brocca et al., 2014; Koster et al., 2009). Considering all these limitations, the improvements achieved with the mHM-RE module are encouraging and highlight its potential for mesoscale hydrological modeling.

Figure 12 shows time series of soil moisture ($\theta$) and soil evapotranspiration ($E_i$) at the Aachen validation site. All model variants display pronounced SM fluctuations in the shallow layers (L1–L3), primarily driven by high $E_i$ rates, which can reach up to 4 mm/d in these layers. As $E_i$ decreases with depth, the behavior of mHM-IC diverges from that of the mHM-RE models. In mHM-IC, soil moisture becomes increasingly stable with depth and eventually reaches saturation in the deepest layer, evident from the near-constant $\theta$ time series.

In contrast, the mHM-RE models continue to exhibit SM fluctuations, even in the deeper layers, despite minimal $E_i$ fluxes. For example, mHM-RE 2 (green) and mHM-RE 3 (red) show negligible $E_i$ starting at layer L4, yet maintain dynamic SM behavior down to layer L6. This sustained variability results from the RE implementation, which permits bidirectional soil moisture flow. When the upper layers experience moisture depletion due to $E_i$ losses, upward fluxes from deeper layers are generated, enabling continued SM redistribution throughout the soil profile. This behavior provides a foundation for future work in which mHM can be coupled with a fully distributed groundwater model.

# 4. Summary and Conclusions

This study evaluated three pedo-transfer functions (PTFs) to parameterize the one-dimensional RE (1D RE) within the mHM mHM (mHM-RE), using the MPR technique. Model performance was benchmarked against the original mHM configuration employing the IC approach (mHM-IC), using both streamflow and soil moisture
(SM) observations. Across 201 basins, all models showed comparable streamflow performance based on KGE scores, demonstrating that the integration of the RE formulation into mHM does not compromise streamflow simulation skill.

Soil moisture validation reveals key differences between the RE and IC approaches. All model variants were able to reproduce temporal SM anomalies effectively. However, mHM-IC underperformed in deeper soil layers due to its restriction to downward-only water movement, which led to less dynamic SM responses. For absolute SM, all models showed systematic underestimation and overestimation, respectively—likely a result of calibrating solely against streamflow observations. Among all RE variants, mHM-RE 3 showed the best overall performance, balancing SM dynamics and realism in absolute values most effectively.

Overall, this study demonstrates that by employing an appropriate parameterization—specifically, the MPR technique—the RE approach can yield transferable parameters and robust simulations of both streamflow and soil moisture dynamics in global-scale modeling. This finding implies that LSMs with poor streamflow performance, despite employing RE-based soil moisture dynamics, likely suffer from inadequate parameterization of the RE, a lack of explicit runoff generation processes, or both. A primary limitation of this study is the exclusive use of streamflow data and a single-basin procedure during the calibration, which may result in suboptimal performance for simulating absolute soil moisture and reduce the effectiveness of cross-basin validation. Future research should focus on multivariate calibration strategies that integrates spatially distributed observations, such as evapotranspiration, long-term groundwater levels, soil moisture, and total water storage to more effectively constrain the infiltration-related parameters as well as the runoff-generation processes.

## Data Availability Statement

SRTM DEM data set (Jarvis et al., 2008) can be downloaded from https://srtm.csi.cgiar.org/srtmdata/. CORINE land cover data (EEA, 2020) are accessible via the Copernicus portal at https://land.copernicus.eu/en/products/corine-land-cover. SoilGrids data set (Poggio et al., 2021) is available at https://soilgrids.org/. MODIS LAI data set (Myneni et al., 2015) can be accessed via http://doi.org/10.5067/MODIS/MYD15A2H.006. GLIM lithology data set (Hartmann & Moosdorf, 2012a, 2012b) is downloadable at https://doi.pangaea.de/10.1594/PANGAEA.788537. German Weather Service (DWD) meteorological data are provided through the DWD portal at https://dwd-geoportal.de/. Streamflow observations from the Global Runoff Data Centre (GRDC) are available on the GRDC data portal at https://portal.grdc.bafg.de/. The original soil moisture (SM) observations from TERENO (Zacharias et al., 2024) and FLUXNET (Pastorello et al., 2020) are available at https://ddp.tereno.net/ddp/ and https://fluxnet.org/data/fluxnet2015-dataset/, respectively. The CRNS SM measurement (Bogena et al., 2022; Bogena & Ney, 2021) can be obtained via Forschungszentrum Jülich data portal (https://doi.org/10.34731/x9s3-kr48). The depth averaged SM data from SPM, SDM and LYSI measurements (Boeing, 2024) are accessible at https://doi.org/10.48758/ufz.15298. The mHM source code for simulating mHM-IC and all mHM-RE variants is available at https://git.ufz.de/afid.nur-kholis/mhm_sli (open access), as a git fork of the main mHM repository (Samaniego et al., 2023). The Python scripts used for generating plots and maps are accessible at https://git.ufz.de/afid.nur-kholis/mhm_re_paper (open access).

## References

Beck, H. E., van Dijk, A. I. J. M., de Roo, A., Miralles, D. G., McVicar, T. R., Schellekens, J., & Bruijnzeel, L. A. (2016). Global-scale regionalization of hydrologic model parameters. *Water Resources Research*, 52(5), 3599–3622. https://doi.org/10.1002/2015WR018247

Bergström, S. (1995). The HBV model. In V. Singh (Ed.), *Computer models of watershed hydrology* (pp. 443–476). Highlands Ranch, Colorado: Water Resources Publications.

Beven, K. (1989). Changing ideas in hydrology—The case of physically-based models. *Journal of Hydrology*, 105(1), 157–172. https://doi.org/10.1016/0022-1694(89)90101-7

Bittelli, M., Campbell, G. S., & Tomei, F. (2015). *Soil physics with python: Transport in the soil-plant-atmosphere system* (1st ed.). Oxford University Press.

Blöschl, G. B., Marc, F. P. C., Antonio, C., Christophe, D., Georgia, F., Fiori, A., et al. (2019). Twenty-three unsolved problems in hydrology (UPH)—A community perspective. *Hydrological Sciences Journal*, 64(10), 1141–1158. https://doi.org/10.1080/02626667.2019.1620507

Boeing, F. (2024). Depth-averaged soil moisture observations at 26 locations in Germany. *Helmholtz-Zentrum für Umweltforschung*. https://doi.org/10.48758/UFZ.15298

Boeing, F., Rakovec, O., Kumar, R., Samaniego, L., Schrön, M., Hildebrandt, A., et al. (2022). High-resolution drought simulations and comparison to soil moisture observations in Germany. *Hydrology and Earth System Sciences*, 26(19), 5137–5161. https://doi.org/10.5194/hess-26-5137-2022

Bogena, H., & Ney, P. (2021). Dataset of COSMOS-Europe: A European network of cosmic-ray neutron soil moisture sensors. *Forschungszentrum Jülich*. https://doi.org/10.34731/x9s3-kr48

## Acknowledgments

This work is part of the Helmholtz Centre for Environmental Research—UFZ internal project LISA (Landtrans Interface and Scaling Solutions, https://www.ufz.de/index.php?en=49789). We would like to thank Pallav K. Shrestha for valuable discussions regarding the mHM model setup. We would like to thank the administration and support staff of EVE (High-Performance Computing Cluster, a joint effort of both UFZ and the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig) who keep the system running and support us with our scientific computing needs: Thomas Schnicke, Ben Langenberg, Guido Schramm, Toni Harzendorf, Tom Strempel, Lisa Schurack, Mark Fliak and Conrad Ostertag from the UFZ, and Christian Krause from iDiv. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). We acknowledge the TERENO community and FLUXNET community for providing soil moisture observational data. We kindly acknowledge the German Weather Service (DWD), the European Environmental Agency (EEA), the European Space Agency (ESA), the U.S. Geological Survey (USGS), International Soil Reference and Information Centre (ISRIC—World Soil Information), and the Global Runoff Data Centre (GRDC) as data providers. Open Access funding enabled and organized by Projekt DEAL.
Bogena, H., Schrön, M., Jakobi, J., Ney, P., Zacharias, S., Andreasen, M., et al. (2022). COSMOS-Europe: A European network of cosmic-ray neutron soil moisture sensors. *Earth System Science Data*, 14(3), 1125-1151. https://doi.org/10.5194/essd-14-1125-2022

Brocca, L., Zucco, G., Mittelbach, H., Moramarco, T., & Seneviratne, S. I. (2014). Absolute versus temporal anomaly and percent of saturation soil moisture spatial variability for six networks worldwide. *Water Resources Research*, 50(7), 5560–5576. https://doi.org/10.1002/2014WR015684

Brooks, R. H., & Corey, A. T. (1964). Hydraulic properties of porous media. *Hydrology Papers*, 3, 1–27.

Campbell, G. S. (1974). A simple method for determining unsaturated conductivity from moisture retention data. *Soil Science*, 117(6), 311–314. https://doi.org/10.1097/00010694-197406000-00001

Campbell, G. S. (1985). *Soil physics with basic: Transport models for soil-plant systems (No. 14)*. Elsevier.

Campbell, G. S., & Shiozawa, S. (1992). Prediction of hydraulic properties of soils using particle size distribution and bulk density data. In *International workshop on indirect methods for estimating the hydraulic properties of unsaturated soils* (pp. 317–328). University of California.

Clapp, R. B., & Hornberger, G. M. (1978). Empirical equations for some soil hydraulic properties. *Water Resources Research*, 14(4), 601–604. https://doi.org/10.1029/WR014i004p00601

Clark, M. P., Fan, Y., Lawrence, D. M., Adam, J. C., Bolster, D., Gochis, D. J., et al. (2015). Improving the representation of hydrologic processes in Earth system models. *Water Resources Research*, 51(8), 5929–5956. https://doi.org/10.1002/2015WR017096

Cosby, B. J., Hornberger, G. M., Clapp, R. B., & Ginn, T. R. (1984). A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils. *Water Resources Research*, 20(6), 682–690. https://doi.org/10.1029/WR020i006p00682

Cuntz, M., & Haverd, V. (2018). Physically accurate soil freeze-thaw processes in a global land surface scheme. *Journal of Advances in Modeling Earth Systems*, 10(1), 54–77. https://doi.org/10.1002/2017MS001100

Cuntz, M., Mai, J., Samaniego, L., Clark, M., Wulfmeyer, V., Branch, O., et al. (2016). The impact of standard and hard-coded parameters on the hydrologic fluxes in the Noah-MP land surface model. *Journal of Geophysical Research: Atmospheres*, 121(18), 10676–10700. https://doi.org/10.1002/2016jd025097

Dai, Y., Xin, Q., Wei, N., Zhang, Y., Shangguan, W., Yuan, H., et al. (2019). A global high-resolution data set of soil hydraulic and thermal properties for land surface modeling. *Journal of Advances in Modeling Earth Systems*, 11(9), 2996–3023. https://doi.org/10.1029/2019MS001784

Dolich, A., Espinoza, E. A., Ebeling, P., Guse, B., Götte, J., Hassler, S., et al. (2024). CAMELS-DE: Hydrometeorological time series and attributes for 1582 catchments in Germany. https://doi.org/10.5281/zenodo.13837553

Döll, P., Douville, H., Güntner, A., Schmied, H. M., & Wada, Y. (2016). Modelling freshwater resources at the global scale: Challenges and prospects. *Surveys in Geophysics*, 37(2), 195–221. https://doi.org/10.1007/s10712-015-9343-1

EEA. (2020). *Corine land cover 1990, 2000, 2006, 2012, 2018 vector/raster 100 m*. Copernicus Land Monitoring Service. Retrieved from https://land.copernicus.eu/en/products/corine-land-cover

Farthing, M. W., & Ogden, F. L. (2017). Numerical solution of Richards' equation: A review of advances and challenges. *Soil Science Society of America Journal*, 81(6), 1257–1269. https://doi.org/10.2136/ssaj2017.02.0058

Feddes, R. A., Kowalik, P., Kolinska-Malinka, K., & Zaradny, H. (1976). Simulation of field water uptake by plants using a soil water dependent root extraction function. *Journal of Hydrology*, 31(1-2), 13–26. https://doi.org/10.1016/0022-1694(76)90017-2

Feigl, M., Herrnegger, M., Klotz, D., & Schulz, K. (2020). Function space optimization: A symbolic regression method for estimating parameter transfer functions for hydrological models. *Water Resources Research*, 56(10), e2028WR027385. https://doi.org/10.1029/2020WR027385

Feigl, M., Thober, S., Schweppe, R., Herrnegger, M., Samaniego, L., & Schulz, K. (2022). Automatic regionalization of model parameters for hydrological models. *Water Resources Research*, 58(12), e2022WR031966. https://doi.org/10.1029/2022WR031966

Gardner, W. R. (1958). Some steady-state solutions to the unsaturated flow equation with application to evaporation from a water table. *Soil Science*, 85(85), 228–232. https://doi.org/10.1097/00010694-195804000-00006

Gleeson, T., Wagener, T., Döll, P., Zipper, S. C., West, C., Wada, Y., et al. (2021). GMD perspective: The quest to improve the evaluation of groundwater representation in continental-to global-scale models. *Geoscientific Model Development*, 14(12), 7545–7571. https://doi.org/10.5194/gmd-14-7545-2021

Göhler, M., Mai, J., & Cuntz, M. (2013). Use of Eigen decomposition in a parameter sensitivity analysis of the community land model. *Journal of Geophysical Research: Biogeosciences*, 118(2), 904–921. https://doi.org/10.1002/jgrg.20072

Guber, A. K., Pachepsky, Y. A., van Genuchten, M. T., Rawls, W. J., Simunek, J., Jacques, D., et al. (2006). Field-scale water flow simulations using ensembles of pedotransfer functions for soil water retention. *Vadose Zone Journal*, 5(1), 234–247. https://doi.org/10.2136/vzj2005.0111

Gupta, H. V., Kling, H., Yilmaz, K. K., & Martinez, G. F. (2009). Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. *Journal of Hydrology*, 377(1), 80–91. https://doi.org/10.1016/j.jhydrol.2009.08.003

Hargreaves, G. H., & Samani, Z. A. (1985). Reference crop evapotranspiration from temperature. *Applied Engineering in Agriculture*, 1(2), 96–99. https://doi.org/10.13031/2013.26773

Hartmann, J., & Moosdorf, N. (2012a). *Global lithological map database v1.0 (gridded to 0.5° spatial resolution)*. PANGAEA. https://doi.org/10.1594/PANGAEA.788537

Hartmann, J., & Moosdorf, N. (2012b). The new global lithological map database GLiM: A representation of rock properties at the Earth surface. *Geochemistry, Geophysics, Geosystems*, 13(12), Q12004. https://doi.org/10.1029/2012GC004370

Haughton, N., Abramowitz, G., Pitman, A. J., Or, D., Best, M. J., Johnson, H. R., et al. (2016). The plumbing of land surface models: Is poor performance a result of methodology or data quality? *Journal of Hydrometeorology*, 17(6), 1705–1723. https://doi.org/10.1175/jhm-d-15-0171.1

Haverd, V., & Cuntz, M. (2010). Soil–litter–ISO: A one-dimensional model for coupled transport of heat, water and stable isotopes in soil with a litter layer and root extraction. *Journal of Hydrology*, 388(3), 438–455. https://doi.org/10.1016/j.jhydrol.2010.05.029

Haverd, V., Cuntz, M., Nieradzik, L. P., & Harman, I. N. (2016). Improved representations of coupled soil–canopy processes in the CABLE land surface model. *Geoscientific Model Development*, 9(9), 3111–3122. https://doi.org/10.5194/gmd-9-3111-2016

Haverd, V., Raupach, M. R., Briggs, P. R., Canadell, J. G., Isaac, P., Pickett-Heaps, C., et al. (2013). Multiple observation types reduce uncertainty in Australia's terrestrial carbon and water cycles. *Biogeosciences*, 10(3), 2011–2040. https://doi.org/10.5194/bg-10-2011-2013

Henderson-Sellers, A., Pitman, A. J., Love, P. K., Irannejad, P., & Chen, T. H. (1995). The project for intercomparison of land surface parameterization schemes (PILPS): Phases 2 and 3. *Bulletin of the American Meteorological Society*, 76(4), 489–504. https://doi.org/10.1175/1520-0477(1995)076%3C0489:TPFIOL%3E2.0.CO;2

Henderson-Sellers, A., Yang, Z.-L., & Dickinson, R. E. (1993). The project for intercomparison of land-surface parameterization schemes. *Bulletin of the American Meteorological Society*, 74(7), 1335–1350. https://doi.org/10.1175/1520-0477(1993)074(1335:TPFIOL)2.0.CO;2
Horton, R. E. (1941). An approach toward a physical interpretation of infiltration-capacity. *Soil Science Society of America Journal*, 5(C), 399-417. https://doi.org/10.2136/sssaj1941.036159950005000C0075x

Hou, Y., Guo, H., Yang, Y., & Liu, W. (2023). Global evaluation of runoff simulation from climate, hydrological and land surface models. *Water Resources Research*, 59(1), e2021WR031817. https://doi.org/10.1029/2021WR031817

Jarvis, A., Reuter, H. I., Nelson, A., & Guevara, E. (2008). Hole-filled seamless SRTM data V4. *International Centre for Tropical Agriculture (CIAT)*. Retrieved from https://srtm.csi.cgiar.org

Kirchner, J. W. (2006). Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology. *Water Resources Research*, 42(3), W03S04. https://doi.org/10.1029/2005WR004362

Klemeš, V. (1983). Conceptualization and scale in hydrology. *Journal of Hydrology*, 65(1-3), 1-23. https://doi.org/10.1016/0022-1694(83)90208-1

Klotz, D., Herrnegger, M., & Schulz, K. (2017). Symbolic regression for the estimation of transfer functions of hydrological models. *Water Resources Research*, 53(11), 9402-9423. https://doi.org/10.1002/2017WR021253

Koirala, S., Yeh, P. J.-F., Hirabayashi, Y., Kanae, S., & Oki, T. (2014). Global-scale land surface hydrologic modeling with the representation of water table dynamics. *Journal of Geophysical Research: Atmospheres*, 119(1), 75-89. https://doi.org/10.1002/2013JD020398

Koster, R. D., Dirmeyer, P. A., Guo, Z., Bonan, G., Chan, E., Cox, P., et al. (2004). Regions of strong coupling between soil moisture and precipitation. *Science*, 305(5687), 1138-1140. https://doi.org/10.1126/science.1100217

Koster, R. D., Guo, Z., Yang, R., Dirmeyer, P. A., Mitchell, K., & Puma, M. J. (2009). On the nature of soil moisture in land surface models. *Journal of Climate*, 22(16), 4322-4335. https://doi.org/10.1175/2009JCLI2832.1

Krysanova, V., Zaherpour, J., Didovets, I., Gosling, S. N., Gerten, D., Hanasaki, N., et al. (2020). How evaluation of global hydrological models can help to improve credibility of river discharge projections under climate change. *Climatic Change*, 163(3), 1353-1377. https://doi.org/10.1007/s10584-020-02840-0

Kumar, R. (2010). Distributed hydrologic model parameterization: Application in a mesoscale river basin. *Doctoral Dissertation Friedrich-Schiller-Universitat Jena*. Retrieved from https://www.ufz.de/index.php?en=20939&pub_id=11006

Kumar, R., Livneh, B., & Samaniego, L. (2013a). Toward computationally efficient large-scale hydrologic predictions with a multiscale regionalization scheme. *Water Resources Research*, 49(9), 5700-5714. https://doi.org/10.1002/wcr.20431

Kumar, R., Samaniego, L., & Attinger, S. (2013b). Implications of distributed hydrologic model parameterization on water fluxes at multiple scales and locations. *Water Resources Research*, 49(1), 360-379. https://doi.org/10.1029/2012WR012195

Lehmann, P., Assouline, S., & Or, D. (2008). Characteristic lengths affecting evaporative drying of porous media. *Physical Review E*, 77(5), 056309. https://doi.org/10.1103/PhysRevE.77.056309

Lehmann, P., & Or, D. (2009). Evaporation and capillary coupling across vertical textural contrasts in porous media. *Physical Review E*, 80(4), 046318. https://doi.org/10.1103/PhysRevE.80.046318

Liang, X., Lettenmaier, D. P., Wood, E. F., & Burges, S. J. (1994). A simple hydrologically based model of land surface water and energy fluxes for general circulation models. *Journal of Geophysical Research*, 99(D7), 14415-14428. https://doi.org/10.1029/94JD00483

Liang, X., & Xie, Z. (2003). Important factors in land-atmosphere interactions: Surface runoff generations and interactions between surface and groundwater. *Global and Planetary Change*, 38(1-2), 101-114. https://doi.org/10.1016/S0921-8181(03)00012-2

Mendoza, P. A., Clark, M. P., Barlage, M., Rajagopalan, B., Samaniego, L., Abramowitz, G., & Gupta, H. (2015). Are we unnecessarily constraining the agility of complex process-based models? *Water Resources Research*, 51(1), 716-728. https://doi.org/10.1002/2014wr015820

Miguez-Macho, G., Fan, Y., Weaver, C. P., Walko, R., & Robock, A. (2007). Incorporating water table dynamics in climate modeling: 2. Formulation, validation, and soil moisture simulation. *Journal of Geophysical Research*, 112(D13), D13108. https://doi.org/10.1029/2006JD008112

Mizukami, N., Clark, M. P., Newman, A. J., Wood, A. W., Gutmann, E. D., Nijssen, B., et al. (2017). Towards seamless large-domain parameter estimation for hydrologic models. *Water Resources Research*, 53(9), 8020-8040. https://doi.org/10.1002/2017WR020401

Mizukami, N., Rakovec, O., Newman, A. J., Clark, M. P., Wood, A. W., Gupta, H. V., & Kumar, R. (2019). On the choice of calibration metrics for "high-flow" estimation using hydrologic models. *Hydrology and Earth System Sciences*, 23(6), 2601-2614. https://doi.org/10.5194/hess-23-2601-2019

Myneni, R., Knyazikhin, Y., & Park, T. (2015). *Myd15a2h modis/aqua leaf area index/fpar 8-day 14 global 500m sin grid. NASA lp daac*. Boston University and MODAPS SIPS-NASA. https://doi.org/10.5067/MODIS/MYD15A2H.006

Nearing, G. S., Mocko, D. M., Peters-Lidard, C. D., Kumar, S. V., & Xia, Y. (2016). Benchmarking NLDAS-2 soil moisture and evapotranspiration to separate uncertainty contributions. *Journal of Hydrometeorology*, 17(3), 745-759. https://doi.org/10.1175/jhm-d-15-0063.1

Newman, A. J., Mizukami, N., Clark, M. P., Wood, A. W., Nijssen, B., & Nearing, G. (2017). Benchmarking of a physically based hydrologic model. *Journal of Hydrometeorology*, 18(8), 2215-2225. https://doi.org/10.1175/JHM-D-16-0284.1

Or, D., Lehmann, P., & Assouline, S. (2015). Natural length scales define the range of applicability of the Richards equation for capillary flows. *Water Resources Research*, 51(9), 7130-7144. https://doi.org/10.1002/2015WR017034

Paniconi, C., Aldama, A. A., & Wood, E. F. (1991). Numerical evaluation of iterative and noniterative methods for the solution of the nonlinear Richards equation. *Water Resources Research*, 27(6), 1147-1163. https://doi.org/10.1029/91WR00334

Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., et al. (2020). The FLUXNET2015 dataset and the one flux processing pipeline for eddy covariance data. *Scientific Data*, 7(1), 225. https://doi.org/10.1038/s41597-020-0534-3

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., & Rossiter, D. (2021). SoilGrids 2.0: Producing soil information for the globe with quantified spatial uncertainty. *SOIL*, 7(1), 217-240. https://doi.org/10.5194/soil-7-217-2021

Pringle, M. J., Romano, N., Minasny, B., Chirico, G. B., & Lark, R. M. (2007). Spatial evaluation of pedotransfer functions using wavelet analysis. *Journal of Hydrology*, 333(2), 182-198. https://doi.org/10.1016/j.jhydrol.2006.08.007

Rakovec, O., Kumar, R., Attinger, S., & Samaniego, L. (2016a). Improving the realism of hydrologic model functioning through multivariate parameter estimation. *Water Resources Research*, 52(10), 7779-7792. https://doi.org/10.1002/2016WR019430

Rakovec, O., Kumar, R., Mai, J., Cuntz, M., Thober, S., Zink, M., et al. (2016b). Multiscale and multivariate evaluation of water fluxes and states over European River basins. *Journal of Hydrometeorology*, 17(1), 287-307. https://doi.org/10.1175/JHM-D-15-0054.1

Rakovec, O., Mizukami, N., Kumar, R., Newman, A. J., Thober, S., Wood, A. W., et al. (2019). Diagnostic evaluation of large-domain hydrologic models calibrated across the contiguous United States. *Journal of Geophysical Research: Atmospheres*, 124(24), 13991-14007. https://doi.org/10.1029/2019JD030767

Rawls, W. J., & Brakensiek, D. L. (1989). Estimation of soil water retention and hydraulic properties. In H. J. Morel-Seytoux (Ed.), *Unsaturated flow in hydrologic modeling* (pp. 275-300). Springer Netherlands. https://doi.org/10.1007/978-94-009-2352-2_10

Ross, P. J. (2003). Modeling soil water and solute transport—fast, simplified numerical solutions. *Agronomy Journal*, 95(6), 1352-1361. https://doi.org/10.2134/agronj2003.1352
Samaniego, L. (2025). Permanent shifts in the global water cycle. *Science*, **387**(6741), 1348–1350. https://doi.org/10.1126/science.adw5851

Samaniego, L., & Bárdossy, A. (2005). Robust parametric models of runoff characteristics at the mesoscale. *Journal of Hydrology*, **303**(1–4), 136–151. https://doi.org/10.1016/j.jhydrol.2004.08.022

Samaniego, L., Kumar, R., & Attinger, S. (2010). Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale. *Water Resources Research*, **46**(5), W05523. https://doi.org/10.1029/2008WR007327

Samaniego, L., Kumar, R., & Jackisch, C. (2011). Predictions in a data-sparse region using a regionalized grid-based hydrologic model driven by remotely sensed data. *Hydrology Research*, **42**(5), 338–355. https://doi.org/10.2166/nh.2011.156

Samaniego, L., Kumar, R., Thober, S., Rakovec, O., Zink, M., Wanders, N., et al. (2017). Toward seamless hydrologic predictions across spatial scales. *Hydrology and Earth System Sciences*, **21**(9), 4323–4346. https://doi.org/10.5194/hess-21-4323-2017

Samaniego, L., Kumar, R., Zink, M., Cuntz, M., Mai, J., Thober, S., et al. (2023). mhm-ufz/mhm: V5.13.1. Zenodo. Zenodo, https://doi.org/10.5281/zenodo.8279545

Schweppe, R., Thober, S., Müller, S., Kelbling, M., Kumar, R., Attinger, S., & Samaniego, L. (2022). Mpr 1.0: A stand-alone multiscale parameter regionalization tool for improved parameter estimation of land surface models. *Geoscientific Model Development*, **15**(2), 859–882. https://doi.org/10.5194/gmd-15-859-2022

Seneviratne, S. I., Corti, T., Davin, E. L., Hirschi, M., Jaeger, E. B., Lehner, I., et al. (2010). Investigating soil moisture-climate interactions in a changing climate: A review. *Earth-Science Reviews*, **99**(3–4), 125–161. https://doi.org/10.1016/j.earscirev.2010.02.004

Shrestha, P. K., Samaniego, L., Rakovec, O., Kumar, R., Mi, C., Rinke, K., & Thober, S. (2024). Toward improved simulations of disruptive reservoirs in global hydrological modeling. *Water Resources Research*, **60**(4), e2023WR035433. https://doi.org/10.1029/2023WR035433

Smith, M. B., Koren, V., Reed, S., Zhang, Z., Zhang, Y., Moreda, F., et al. (2012). The distributed model intercomparison project-Phase 2: Motivation and design of the Oklahoma experiments. *Journal of Hydrology*, **418–419**, 3–16. https://doi.org/10.1016/j.jhydrol.2011.08.055

Smith, M. B., Seo, D.-J., Koren, V. I., Reed, S. M., Zhang, Z., Duan, Q., et al. (2004). The distributed model intercomparison project (DMIP): Motivation and experiment design. *Journal of Hydrology*, **298**(1), 4–26. https://doi.org/10.1016/j.jhydrol.2004.03.040

Soylu, M. E., Istambulluoglu, E., Lenters, J. D., & Wang, T. (2011). Quantifying the impact of groundwater depth on evapotranspiration in a semi-arid grassland region. *Hydrology and Earth System Sciences*, **15**(3), 787–806. https://doi.org/10.5194/hess-15-787-2011

Thober, S., Cuntz, M., Kelbling, M., Kumar, R., Mai, J., & Samaniego, L. (2019). The multiscale routing model mRM v1.0: Simple river routing at resolutions from 1 to 50 km. *Geoscientific Model Development*, **12**(6), 2501–2521. https://doi.org/10.5194/gmd-12-2501-2019

Tolson, B. A., & Shoemaker, C. A. (2007). Dynamically dimensioned search algorithm for computationally efficient watershed model calibration. *Water Resources Research*, **43**(1), W01413. https://doi.org/10.1029/2005WR004723

Twarakavi, N. K. C., Sakai, M., & Šimáunek, J. (2009). An objective analysis of the dynamic nature of field capacity. *Water Resources Research*, **45**(10), W10410. https://doi.org/10.1029/2009WR007944

van Beek, L. P. H., & Bierkens, M. F. P. (2009). The global hydrological model PCR-GLOBWB: Conceptualization, parameterization and verification (Technical Reports). *Department of Physical Geography*, Retrieved from https://vanbeek.geo.uu.nl/suppinfo/vanbeekbierkens2009.pdf

van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Science Society of America Journal*, **44**(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x

van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., et al. (2017). Pedotransfer functions in Earth system science: Challenges and perspectives. *Reviews of Geophysics*, **55**(4), 1199–1256. https://doi.org/10.1002/2017RG000581

Vereecken, H., Weihermüller, L., Assouline, S., Šimáunek, J., Verhoef, A., Herbst, M., et al. (2019). Infiltration from the pedon to global grid scales: An overview and outlook for land surface modeling. *Vadose Zone Journal*, **18**(1), 1–53. https://doi.org/10.2136/vzj2018.10.0191

Vogel, H.-J., & Ippisch, O. (2008). Estimation of a critical spatial discretization limit for solving Richards’ equation at large scales. *Vadose Zone Journal*, **7**(1), 112–114. https://doi.org/10.2136/vzj2006.0182

Warszawski, L., Frieler, K., Huber, V., Piontek, F., Serdeczny, O., & Schewe, J. (2014). The inter-sectoral impact model intercomparison project (ISI-MIP): Project framework. *Proceedings of the National Academy of Sciences*, **111**(9), 3228–3232. https://doi.org/10.1073/pnas.1312330110

Weihermüller, L., Lehmann, P., Herbst, M., Rahmati, M., Verhoef, A., Or, D., et al. (2021). Choice of pedotransfer functions matters when simulating soil water balance fluxes. *Journal of Advances in Modeling Earth Systems*, **13**(3), e2020MS002404. https://doi.org/10.1029/2020MS002404

Wood, E. F., Roundy, J. K., Troy, T. J., van Beek, L. P. H., Bierkens, M. F. P., Blyth, E., et al. (2011). Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water. *Water Resources Research*, **47**(5), W05301. https://doi.org/10.1029/2010WR010090

Zacharias, S., Loescher, H. W., Bogena, H., Kiese, R., Schrön, M., Attinger, S., et al. (2024). Fifteen years of integrated terrestrial environmental observatories (TERENO) in Germany: Functions, services, and lessons learned. *Earth’s Future*, **12**(6), e2024EF004510. https://doi.org/10.1029/2024EF004510

Zacharias, S., & Wessolek, G. (2007). Excluding organic matter content from pedotransfer predictors of soil water retention. *Soil Science Society of America Journal*, **71**(1), 43–50. https://doi.org/10.2136/sssaj2006.0098

Zhang, Y., Schaap, M. G., & Wei, Z. (2020). Development of hierarchical ensemble model and estimates of soil water retention with global coverage. *Geophysical Research Letters*, **47**(15), e2020GL088819. https://doi.org/10.1029/2020GL088819

Zink, M., Kumar, R., Cuntz, M., & Samaniego, L. (2017). A high-resolution dataset of water fluxes and states for Germany accounting for parametric uncertainty. *Hydrology and Earth System Sciences*, **21**(3), 1769–1790. https://doi.org/10.5194/hess-21-1769-2017

Zink, M., Mai, J., Cuntz, M., & Samaniego, L. (2018). Conditioning a hydrologic model using patterns of remotely sensed land surface temperature. *Water Resources Research*, **54**(4), 2976–2998. https://doi.org/10.1002/2017WR021346
# Supporting Information for "Evaluating Richards equation and infiltration capacity approaches in mesoscale hydrologic modelling"

A. Kholis<fnref n="1" />,<fnref n="3" />, T. Kalbacher<fnref n="1" />, O. Rakovec<fnref n="2" />,<fnref n="4" />, F. Boeing<fnref n="2" />,<fnref n="3" />, M. Cuntz<fnref n="5" />, L.

Samaniego<fnref n="2" />,<fnref n="3" />

## Contents of this file

1. Text S1 to S5

2. Figures S1 to S11

3. Tables S1 to S10

## Introduction

The following supporting information provides a description of the main data used for model setup and validation (Text S1). An overview of the global parameters (γ) involved in the calibration process, including their ranges and initial values, is detailed in Text S2. Additionally, supplementary figures for the cross-basin validation (Text S3),

<footnote n="1">¹Department Environmental Informatics, UFZ - Helmholtz Centre for Environmental Research GmbH, Leipzig, Germany</footnote>

<footnote n="2">²Department Computational Hydrosystems, UFZ - Helmholtz Centre for Environmental Research GmbH, Leipzig, Germany</footnote>

<footnote n="3">³University of Potsdam, Institute of Environmental Science and Geography, Potsdam, Germany</footnote>

<footnote n="4">⁴Faculty of Environmental Sciences, Czech University of Life Sciences, Prague, Czech Republic</footnote>

<footnote n="5">⁵Université de Lorraine, AgroParisTech, INRAE, UMR Silva, Nancy, France</footnote>
Germany-wide simulations (Text S4), and soil moisture validation (Text S5) are included
and explained in this text.

**Text S1. Study Area and Datasets**

Details of the datasets used in this study are provided in Table S1. The model's under-
lying setup data was constructed using a base l₀ grid with a 200 m × 200 m resolution.
Elevation data from the SRTM DEM (Jarvis et al., 2008) was pre-processed using GIS
software to derive additional terrain characteristics, including slope, aspect, flow direction,
and flow accumulation. The CORINE dataset (EEA, 2020), which includes land cover
data for five different years (1990, 2000, 2006, 2012, and 2018), was reclassified into three
land cover types supported by the mHM model: forest, permeable areas, and impervious
cover.

The climatology of leaf area index (LAI) was derived from MODIS Terra version 6
(Myneni et al., 2015), comprising mean monthly LAI values for the period 2003–2018.
Soil physical characteristics (sand content, clay content, and bulk density) were obtained
from the SoilGrids dataset (Poggio et al., 2021), available at six depths: 0-5 cm (L1),
5-15 cm (L2), 15-25 cm (L3), 25-60 cm (L4), 60-100 cm (L5), and 100-200 cm (L6).
Lithology data from the GLIM dataset (Hartmann & Moosdorf, 2012a, 2012b) covers 16
different lithological types, and the lithology present in each basin was calibrated as GeoP-
aram parameters (Table S3). The meteorological forcing data (precipitation and mean air
temperature *T*<sub>mean</sub>) available at 1 km × 1 km resolution were derived from the German
Weather Service (DWD) gauges using external drift kriging (EDK) conditioned on a dig-
ital elevation model (Samaniego et al., 2011). The potential evapotranspiration (PET)
was calculated using the Hargreaves and Samani (1985) method, based on maximum and minimum daily air temperature (*T*<sub>max</sub> and *T*<sub>min</sub>) obtained from the DWD.

The reference soil moisture (SM) observation data at plot scale, covering 46 sites (0–25 cm depth) and 42 sites (25–60 cm and 0–60 cm depths) across Germany, was obtained from Boeing et al. (2022) (Table S2). The data originally comes from FLUXNET observation network (Pastorello et al., 2020) and four TERENO sites (Zacharias et al., 2024): Central Germany (CG), Rur/Eifel (Rur/E), Northeast (NE), and the Pre-Alpine Observatory (PAO). Figure 2 in the main text shows the distribution of these measurement networks, which span five land use types: cropland, grassland, deciduous broadleaf forest (DBF), clearings, and evergreen needleleaf forest (ENF). The SM data covers various time periods and is collected using four methods: cosmic-ray neutron sensing (CRNS), spatially distributed measurements (SDM), single profile measurements (SPM), and lysimeter (LYSI). The original dataset included multiple sensors at different depths at each site, which were aggregated using depth-averaging to align with the model’s soil layer configuration: 0–25 cm (shallow soil, L1–L3), 25–60 cm (deep soil, L4), and 0–60 cm (average of shallow and deep soil layers, L1–L4).

**Text S2. Calibrated Parameters (Global Parameters)**

All global parameters (γ) defining the mHM processes were calibrated in this study. Table S3 lists the parameters for each process—except the soil infiltration—involved in the optimization process for both mHM-IC and mHM-REs variants. The range and initial values were obtained from the mHM's official documentation (https://mhm-ufz.org/) and based on extensive sensitivity studies across various basins.
The primary difference between mHM-IC and mHM-REs lies in the global parameters $\gamma$ governing the soil infiltration process. The $\gamma$ parameters of the infiltration process in mHM-IC (Table S4) are used to calculate porosity $\theta_s$ ($\gamma_1-\gamma_6$), pore size distribution index $\lambda$ ($\gamma_7$), and saturated hydraulic conductivity $K_s$ ($\gamma_8-\gamma_{10}$). Meanwhile, for the mHM-RE variants, the SM global parameters $\gamma$ were used to calculate $K_s$, $\theta_s$, $\lambda$, $h_s$ (air entry value), and, in the case of mHM-RE 3, $\theta_r$ (residual water content), as provided in Tables S5-S7.

The initial values for these infiltration global parameters were derived from the original pedo-transfer function (PTF) equations for each model: Cosby (Cosby et al., 1984) for mHM-RE 1 (Table S5), Campbell (1985) for the mHM-RE 2 model (Table S6), and Rawls and Brakensiek (1989) for the mHM-RE 3 model (Table S7). The calibration range values were subjectively assigned based on the initial values to provide broad parameter exploration for the models.

## Text S3. Cross-basin validation

Figure S1 extends Figure 5, showing the cumulative frequency distribution from the cross-basin validation with 201 domains for each model. The individual runs, based on the optimized parameters from individual basin calibration, are shown in gray. The best simulation for each model, determined by the highest median KGE value, is displayed with a colored line: basin 121 (blue color) for mHM-IC, basin 120 (orange color) for mHM-RE 1, basin 171 (green color for mHM-RE 2), and basin 134 (red color) for mHM-RE 3.

## Text S4. Germany-wide simulation

Figures S2 and S3 expand on the mean runoff maps shown in Figure 8. All maps represent the average values over the 30-year simulation period (1989-2018). Figure S2
provides the spatial distributions of runoff components for each model, including $Q_D$ (direct runoff) from sealed area storage, $Q_{If}$ (fast interflow) and $QIs$ (slow interflow) from subsurface storage, and $Q_B$ (baseflow) from groundwater storage. It is important to note that interflow components in the models (mHM-IC and mHM-RE) are not directly generated from the soil layers; instead, they are calculated from the subsurface bucket, which is supplied by vertical flux from the last soil layer ($q_n$). The transmission of flux from unsaturated to saturated layers occurs via groundwater recharge (recharge). The mean values of recharge, subsurface storage, and groundwater storage over the simulation period are shown in Figure S3.

**Text S5. Plot scale soil moisture (SM) simulation**

Comparisons between models and observed absolute volumetric SM at 46 sites (0–25 cm depth) and 42 sites (25–60 cm and 0–60 cm depths) are shown in Figures S6–S8, which supplement Figures 11 and 12 in the main text. Comparisons of modeled and observed SM anomalies ($\Delta\theta$) are provided in Figures S9–S11. These figures cover three depth ranges: 0–25 cm (L1 to L3 in the model setup), 25–60 cm (L4 in the model setup), and 0–60 cm (L1 to L4 in the model setup). The model values for the 0–25 cm and 0–60 cm depths were calculated using the depth-averaging method, as they consist of multiple soil layers. The performance evaluation of absolute SM was conducted using KGE and its components (Figure 11 in the main text). On the other hand, the anomaly performance was evaluated using only the Pearson correlation coefficient ($r$), as shown in Figure S5.
<figure><img src="image_15.png" /></figure>

**Figure S1.** Cross-basin validation across 201 domains for each model. The colored line is the best simulation obtained from the highest median KGE value. The number in the legend (e.g., run-121) shows that the model was run with parameter set from that basin
<figure><img src="image_16.png" /><figcaption>**Figure S2.** Maps of Germany-wide models mean runoff (*Q*) components in mm/month: (a)-(d) direct runoff *QD*, (e)-(h) fast interflow *QIf*, (i)-(l) slow interflow *QIs*, and (m)-(p) baseflow *QB*.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure S3.** Maps of Germany-wide models: (a)-(d) groundwater (GW) recharge, (e)-(h) unsaturated (Unsat.) storage, and (i)-(l) saturated (Sat.) storage.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Figure S4.** Maps of Germany-wide models: (a)-(d) interception, (e)-(h) throughflow, (i)-(l) snowpack QIs, and (m)-(p) effective precipitation $P_{eff}$.

June 26, 2025, 11:27pm
<figure><img src="image_19.png" /></figure>

**Figure S5.** Pearson correlation coefficient *r* on soil moisture anomaly (Δθ) over 46 sites (0-25 cm depth) and 42 sites (25-60 cm and 0-60 cm depths) at a daily time step and 200 m grid resolution.
<figure><img src="image_20.png" /><figcaption>**Figure S6.** Observations versus models of volumetric SM ($\theta$) at 0-25 cm depth over 46 sites with a daily time step.</figcaption></figure>

June 26, 2025, 11:27pm
<figure><img src="image_21.png" /><figcaption>**Figure S7.** Observations versus models of volumetric SM ($\theta$) at 25-60 cm depth over 42 sites with a daily time step.</figcaption></figure>

June 26, 2025, 11:27pm
<figure><img src="image_22.png" /></figure>

**Figure S8.** Observations versus models of volumetric SM ($\theta$) at 0-60 cm depth over 42 sites with a daily time step.

June 26, 2025, 11:27pm
<figure><img src="image_23.png" /><figcaption>**Figure S9.** Observations versus models of anomaly SM ($\Delta\theta$) at 0-25 cm depth over 46 sites with a daily time step.</figcaption></figure>

June 26, 2025, 11:27pm
<figure><img src="image_24.png" /></figure>

**Figure S10.** Observations versus models of anomaly SM ($\Delta\theta$) at 25-60 cm depth over 42 sites with a daily time step.

June 26, 2025, 11:27pm
<figure><img src="image_25.png" /></figure>

**Figure S11.** Observations versus models of anomaly SM ($\Delta\theta$) at 0-60 cm depth over 42 sites with a daily time step.

June 26, 2025, 11:27pm
<table><caption>**Table S1.** Summary of datasets used for model setup and evaluation</caption><thead><tr><th>Variable</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td colspan="3"><strong>Model setup</strong></td></tr><tr><td>Meteorological forcing (precipitation, T<sub>mean</sub>, T<sub>min</sub>, and T<sub>max</sub>)</td><td>Daily DWD product at 1 km × 1 km</td><td>Interpolation of point data following Samaniego et al. (2011)</td></tr><tr><td>Terrain characteristics (elevation, slope, aspect, flow direction, and flow accumulation)</td><td>SRTM DEM at 90 m × 90 m</td><td>Jarvis et al. (2008)</td></tr><tr><td>Land cover (forest, permeable, and impervious cover)</td><td>CORINE land-cover at 100 m × 100 m</td><td>EEA (2020)</td></tr><tr><td>Soil properties (sand content, clay content, bulk density, and horizon depth)</td><td>SoilGrids dataset, at 250 m × 250 m</td><td>Poggio et al. (2021)</td></tr><tr><td>LAI</td><td>MODIS dataset, at 500 m × 500 m</td><td>Myneni et al. (2015)</td></tr><tr><td>Lithology</td><td>GLIM dataset, at 0.5° × 0.5°</td><td>Hartmann and Moosdorf (2012a)</td></tr><tr><td colspan="3"><strong>Model evaluation</strong></td></tr><tr><td>Streamflow</td><td>Daily observations at 201 basins</td><td>? (?)</td></tr><tr><td>SM</td><td>Daily observations at 46 sites (0-25 cm depth), 42 sites (25-60 cm and 0-60 cm depths)</td><td>see Table S2</td></tr></tbody></table>
<table><caption>**Table S2.** Overview of soil moisture (SM) measurement sites, retrieved from Boeing et al. (2022)</caption><thead><tr><th>Network</th><th>Site</th><th>Method</th><th>Begin</th><th>End</th><th>Land use</th></tr></thead><tbody><tr><td rowspan="12">TERENO CG</td><td>Bad Lauchstädt</td><td>LYSI</td><td>01/01/2016</td><td>31/12/2018</td><td>Crop</td></tr><tr><td>Ermsleben</td><td>SPM</td><td>25/01/2012</td><td>31/12/2019</td><td>Grass</td></tr><tr><td>Hecklingen</td><td>SPM</td><td>05/07/2013</td><td>31/12/2019</td><td>Grass</td></tr><tr><td>Am Grossen Bruch 1</td><td>CRNS</td><td>24/06/2014</td><td>28/11/2019</td><td>Grass</td></tr><tr><td>Am Grossen Bruch 2</td><td>SDM</td><td>30/07/2014</td><td>18/11/2019</td><td>Grass</td></tr><tr><td>Am Grossen Bruch 3</td><td>SPM</td><td>07/02/2014</td><td>31/12/2019</td><td>Grass</td></tr><tr><td>Hohes Holz 1</td><td>CRNS</td><td>27/08/2014</td><td>28/11/2019</td><td>DBF</td></tr><tr><td>Hohes Holz 2</td><td>SDM</td><td>20/07/2012</td><td>30/12/2019</td><td>DBF</td></tr><tr><td>Hohes Holz 3</td><td>SPM</td><td>25/04/2013</td><td>31/12/2019</td><td>DBF</td></tr><tr><td>Hordorf 1</td><td>CRNS</td><td>29/09/2016</td><td>28/11/2019</td><td>Crop</td></tr><tr><td>Hordorf 2</td><td>SPM</td><td>06/11/2015</td><td>31/12/2019</td><td>Crop</td></tr><tr><td>Cunnersdorf</td><td>CRNS</td><td>23/06/2016</td><td>31/12/2019</td><td>Crop</td></tr><tr><td rowspan="15">TERENO Rur/E</td><td>Aachen</td><td>CRNS</td><td>13/01/2012</td><td>01/05/2019</td><td>Crop</td></tr><tr><td>Gevenich</td><td>CRNS</td><td>06/07/2011</td><td>04/01/2019</td><td>Crop</td></tr><tr><td>Heinsberg</td><td>CRNS</td><td>08/09/2011</td><td>01/05/2019</td><td>Grass</td></tr><tr><td>Kall</td><td>CRNS</td><td>14/09/2011</td><td>01/05/2019</td><td>Grass</td></tr><tr><td>Kleinau</td><td>CRNS</td><td>25/08/2015</td><td>26/04/2019</td><td>Grass</td></tr><tr><td>Merzenhausen</td><td>CRNS</td><td>18/05/2011</td><td>03/04/2019</td><td>Crop</td></tr><tr><td>Rollesbr1</td><td>CRNS</td><td>18/05/2011</td><td>31/12/2018</td><td>Grass</td></tr><tr><td>Rollesbr2</td><td>CRNS</td><td>30/06/2012</td><td>25/12/2018</td><td>Grass</td></tr><tr><td>Ruraue</td><td>CRNS</td><td>08/11/2011</td><td>01/01/2019</td><td>Grass</td></tr><tr><td>Schoeneseiffen</td><td>CRNS</td><td>13/08/2015</td><td>25/04/2019</td><td>Grass</td></tr><tr><td>Selhausen</td><td>CRNS</td><td>06/03/2015</td><td>26/04/2019</td><td>Crop</td></tr><tr><td>Wildenrath</td><td>CRNS</td><td>11/05/2012</td><td>23/03/2019</td><td>Clear</td></tr><tr><td>Wüstebach 1</td><td>CRNS</td><td>12/03/2011</td><td>05/10/2018</td><td>ENF</td></tr><tr><td>Wüstebach 2</td><td>SDM</td><td>27/01/2009</td><td>31/12/2019</td><td>ENF</td></tr><tr><td rowspan="10">TERENO NE</td><td>Alt Tellin</td><td>SPM</td><td>10/05/2014</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Bentzin</td><td>SPM</td><td>19/08/2013</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Droennewitz</td><td>SPM</td><td>12/04/2014</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Goermin</td><td>SPM</td><td>19/08/2013</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Leppin</td><td>SPM</td><td>28/01/2013</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Medrow</td><td>SPM</td><td>13/07/2015</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Muehlenkamp</td><td>SPM</td><td>01/01/2012</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Ueckeritz</td><td>SPM</td><td>28/01/2013</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Wotenick</td><td>SPM</td><td>30/04/2014</td><td>30/12/2019</td><td>Grass</td></tr><tr><td>Zarnekla</td><td>SPM</td><td>23/01/2013</td><td>30/12/2019</td><td>Grass</td></tr><tr><td rowspan="3">TERENO PAO</td><td>Fendt</td><td>LYSI</td><td>01/01/2017</td><td>31/12/2019</td><td>Grass</td></tr><tr><td>Graswang</td><td>LYSI</td><td>17/03/2017</td><td>31/12/2019</td><td>Grass</td></tr><tr><td>Graswang</td><td>LYSI</td><td>17/03/2017</td><td>31/12/2019</td><td>Grass</td></tr><tr><td rowspan="3">FLUXNET</td><td>Gebesee</td><td>SPM</td><td>16/01/2001</td><td>31/12/2014</td><td>Crop</td></tr><tr><td>Grillenburg*</td><td>SPM</td><td>21/11/2006</td><td>31/12/2014</td><td>Grass</td></tr><tr><td>Hainich</td><td>SPM</td><td>27/12/2002</td><td>31/12/2012</td><td>DBF</td></tr></tbody></table>

June 26, 2025, 11:27pm
<table><caption>**Table S2** continued from previous page</caption><thead><tr><th>Network</th><th>Site</th><th>Method</th><th>Begin</th><th>End</th><th>Land use</th></tr></thead><tbody><tr><td></td><td>Klingenberg*</td><td>SPM</td><td>27/11/2004</td><td>31/12/2014</td><td>Crop</td></tr><tr><td></td><td>Lackenberg*</td><td>SPM</td><td>01/05/2009</td><td>31/12/2013</td><td>DBF</td></tr><tr><td></td><td>Leinefelde</td><td>SPM</td><td>01/05/2002</td><td>31/12/2012</td><td>DBF</td></tr><tr><td></td><td>Tharandt*</td><td>SPM</td><td>06/03/1997</td><td>31/12/2014</td><td>DBF</td></tr></tbody></table>

DBF is deciduous broadleaf forest, ENF is evergreen needleleaf forest, *only available at 0-25 cm depth
<table><caption>**Table S3.** Global calibration parameters (γ) for the original mHM model covering all mHM processes (except the infiltration) optimized in this study</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td rowspan="3">Soil organic matter</td><td>γ<sub>1</sub></td><td>orgMatterContent_forest</td><td>0.00</td><td>20.00</td><td>3.40</td></tr><tr><td>γ<sub>2</sub></td><td>orgMatterContent impervious</td><td>0.00</td><td>1.00</td><td>0.10</td></tr><tr><td>γ<sub>3</sub></td><td>orgMatterContent pervious</td><td>0.00</td><td>4.00</td><td>0.60</td></tr><tr><td rowspan="3">Soil evapotranspiration</td><td>γ<sub>1</sub></td><td>rootFractionCoefficient_forest</td><td>0.900</td><td>0.999</td><td>0.970</td></tr><tr><td>γ<sub>2</sub></td><td>rootFractionCoefficient impervious</td><td>0.900</td><td>0.950</td><td>0.930</td></tr><tr><td>γ<sub>3</sub></td><td>rootFractionCoefficient pervious</td><td>0.001</td><td>0.090</td><td>0.020</td></tr><tr><td>Interception</td><td>γ<sub>1</sub></td><td>canopyInterceptionFactor</td><td>0.15</td><td>0.40</td><td>0.15</td></tr><tr><td rowspan="8">Snow</td><td>γ<sub>1</sub></td><td>snowTreshholdTemperature</td><td>-2.00</td><td>2.00</td><td>1.00</td></tr><tr><td>γ<sub>2</sub></td><td>degreeDayFactor_forest</td><td>1.0E-4</td><td>4.00</td><td>1.50</td></tr><tr><td>γ<sub>3</sub></td><td>degreeDayFactor impervious</td><td>0.00</td><td>1.00</td><td>0.50</td></tr><tr><td>γ<sub>4</sub></td><td>degreeDayFactor pervious</td><td>0.00</td><td>2.00</td><td>0.50</td></tr><tr><td>γ<sub>5</sub></td><td>increaseDegreeDayFactorByPrecip</td><td>0.10</td><td>0.90</td><td>0.50</td></tr><tr><td>γ<sub>6</sub></td><td>maxDegreeDayFactor_forest</td><td>0.00</td><td>8.00</td><td>3.00</td></tr><tr><td>γ<sub>7</sub></td><td>maxDegreeDayFactor impervious</td><td>0.00</td><td>8.00</td><td>3.50</td></tr><tr><td>γ<sub>8</sub></td><td>maxDegreeDayFactor pervious</td><td>0.00</td><td>8.00</td><td>4.00</td></tr><tr><td>Direct runoff</td><td>γ<sub>1</sub></td><td>imperviousStorageCapacity</td><td>0.00</td><td>5.00</td><td>0.5</td></tr><tr><td rowspan="5">Interflow</td><td>γ<sub>1</sub></td><td>interflowStorageCapacityFactor</td><td>75.0</td><td>200.0</td><td>85.0</td></tr><tr><td>γ<sub>2</sub></td><td>interflowRecession_slope</td><td>0.00</td><td>10.00</td><td>7.00</td></tr><tr><td>γ<sub>3</sub></td><td>fastInterflowRecession_forest</td><td>1.00</td><td>3.00</td><td>1.50</td></tr><tr><td>γ<sub>4</sub></td><td>slowInterflowRecession_Ks</td><td>1.00</td><td>30.00</td><td>15.00</td></tr><tr><td>γ<sub>5</sub></td><td>exponentSlowInterflow</td><td>0.050</td><td>0.300</td><td>0.125</td></tr><tr><td>Percolation</td><td>γ<sub>1</sub></td><td>rechargeCoefficient</td><td>0.00</td><td>50.00</td><td>35.00</td></tr><tr><td rowspan="16">Baseflow</td><td>γ<sub>1</sub></td><td>baseflowLithology1</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>2</sub></td><td>baseflowLithology2</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>3</sub></td><td>baseflowLithology3</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>4</sub></td><td>baseflowLithology4</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>5</sub></td><td>baseflowLithology5</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>6</sub></td><td>baseflowLithology6</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>7</sub></td><td>baseflowLithology7</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>8</sub></td><td>baseflowLithology8</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>9</sub></td><td>baseflowLithology9</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>10</sub></td><td>baseflowLithology10</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>11</sub></td><td>baseflowLithology11</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>12</sub></td><td>baseflowLithology12</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>13</sub></td><td>baseflowLithology13</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>14</sub></td><td>baseflowLithology14</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>15</sub></td><td>baseflowLithology15</td><td>1.00</td><td>1000</td><td>100</td></tr><tr><td>γ<sub>16</sub></td><td>baseflowLithology16</td><td>1.00</td><td>1000</td><td>100</td></tr></tbody></table>
<table><caption>**Table S4. mHM-IC** model's global calibration parameters ($\gamma$) related to soil moisture (SM) infiltration process. Other mHM processes are parameterized similarly to the original mHM model as defined in Table S3.</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td rowspan="10">Infiltration</td><td>γ<sub>1</sub></td><td>PTF-θ<sub>s1</sub></td><td>0.646</td><td>0.951</td><td>0.760</td></tr><tr><td>γ<sub>2</sub></td><td>PTF-θ<sub>s2</sub></td><td>0.1E-3</td><td>2.9E-3</td><td>0.9E-3</td></tr><tr><td>γ<sub>3</sub></td><td>PTF-θ<sub>s3</sub></td><td>-0.373</td><td>-0.187</td><td>-0.264</td></tr><tr><td>γ<sub>4</sub></td><td>PTF-θ<sub>s4</sub></td><td>0.536</td><td>1.123</td><td>0.890</td></tr><tr><td>γ<sub>5</sub></td><td>PTF-θ<sub>s5</sub></td><td>-5.5E-3</td><td>4.9E-3</td><td>-1.0E-3</td></tr><tr><td>γ<sub>6</sub></td><td>PTF-θ<sub>s6</sub></td><td>-0.551</td><td>-0.091</td><td>-0.324</td></tr><tr><td>γ<sub>7</sub></td><td>infiltrationShapeFactor</td><td>1.00</td><td>4.00</td><td>1.75</td></tr><tr><td>γ<sub>8</sub></td><td>PTF-K<sub>S1</sub></td><td>-1.200</td><td>-0.285</td><td>-0.585</td></tr><tr><td>γ<sub>9</sub></td><td>PTF-K<sub>S2</sub></td><td>6.0E-3</td><td>26.0E-3</td><td>12.5E-3</td></tr><tr><td>γ<sub>10</sub></td><td>PTF-K<sub>S3</sub></td><td>3.0E-3</td><td>13.0E-3</td><td>6.3E-3</td></tr></tbody></table>

$\gamma_1-\gamma_6$ are used to calculate $\theta_s$, $\gamma_7$ is used to calculate $\lambda$, and $\gamma_8-\gamma_{10}$ are used to calculate $K_s$.

<table><caption>**Table S5. mHM-RE 1** model's global calibration parameters ($\gamma$) related to soil moisture (SM) infiltration process. Other mHM processes are parameterized similarly to the original mHM model as defined in Table S3.</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td rowspan="10">Infiltration</td><td>γ<sub>1</sub></td><td>PTF-h<sub>S1</sub></td><td>0.7</td><td>3.00</td><td>1.54</td></tr><tr><td>γ<sub>2</sub></td><td>PTF-h<sub>S2</sub></td><td>-0.02</td><td>-4.0E-3</td><td>-9.5E-3</td></tr><tr><td>γ<sub>3</sub></td><td>PTF-h<sub>S3</sub></td><td>3.0E-3</td><td>12.0E-3</td><td>6.3E-3</td></tr><tr><td>γ<sub>4</sub></td><td>PTF-K<sub>S1</sub></td><td>-1.20</td><td>-0.285</td><td>-0.467</td></tr><tr><td>γ<sub>5</sub></td><td>PTF-K<sub>S2</sub></td><td>0.006</td><td>0.026</td><td>0.016</td></tr><tr><td>γ<sub>6</sub></td><td>PTF-K<sub>S3</sub></td><td>0.003</td><td>0.013</td><td>0.007</td></tr><tr><td>γ<sub>7</sub></td><td>PTF-λ<sub>1</sub></td><td>0.50</td><td>6.00</td><td>3.10</td></tr><tr><td>γ<sub>8</sub></td><td>PTF-λ<sub>2</sub></td><td>0.080</td><td>0.400</td><td>0.157</td></tr><tr><td>γ<sub>9</sub></td><td>PTF-λ<sub>3</sub></td><td>-0.006</td><td>-0.001</td><td>-0.003</td></tr><tr><td>γ<sub>10</sub></td><td>PTF-θ<sub>s1</sub></td><td>30.00</td><td>70.00</td><td>50.50</td></tr><tr><td>γ<sub>11</sub></td><td>PTF-θ<sub>s2</sub></td><td>-0.30</td><td>-0.07</td><td>-0.142</td></tr><tr><td>γ<sub>12</sub></td><td>PTF-θ<sub>s3</sub></td><td>-0.08</td><td>-0.01</td><td>-0.037</td></tr></tbody></table>

$\gamma_1-\gamma_3$ are used to calculate $h_s$, $\gamma_4-\gamma_6$ are used to calculate $K_s$, $\gamma_7-\gamma_9$ are used to calculate $\lambda$, and $\gamma_{10}-\gamma_{12}$ are used to calculate $\theta_s$.
<table><caption>**Table S6.** mHM-RE 2 model's global calibration parameters ($\gamma$) related to soil moisture (SM) infiltration process. Other mHM processes are parameterized similarly to the original mHM model as defined in Table S3.</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td rowspan="10">Infiltration</td><td>γ<sub>1</sub></td><td>PTF-h<sub>S1</sub></td><td>-2.50</td><td>-0.10</td><td>-0.80</td></tr><tr><td>γ<sub>2</sub></td><td>PTF-h<sub>S2</sub></td><td>-0.16</td><td>-0.008</td><td>-0.0317</td></tr><tr><td>γ<sub>3</sub></td><td>PTF-h<sub>S3</sub></td><td>-0.28</td><td>-0.01</td><td>-0.0761</td></tr><tr><td>γ<sub>4</sub></td><td>PTF-K<sub>S1</sub></td><td>1.5E-6</td><td>1.5E-4</td><td>1.5E-5</td></tr><tr><td>γ<sub>5</sub></td><td>PTF-K<sub>S2</sub></td><td>-20.00</td><td>-1.00</td><td>7.00</td></tr><tr><td>γ<sub>6</sub></td><td>PTF-K<sub>S3</sub></td><td>-40.00</td><td>-1.00</td><td>-16.70</td></tr><tr><td>γ<sub>7</sub></td><td>PTF-λ<sub>1</sub></td><td>0.05</td><td>0.50</td><td>0.133</td></tr><tr><td>γ<sub>8</sub></td><td>PTF-λ<sub>2</sub></td><td>0.05</td><td>1.40</td><td>0.477</td></tr><tr><td>γ<sub>9</sub></td><td>PTF-θ<sub>s1</sub></td><td>0.70</td><td>0.999</td><td>0.930</td></tr><tr><td>γ<sub>10</sub></td><td>PTF-θ<sub>s2</sub></td><td>1.00</td><td>5.00</td><td>2.65</td></tr></tbody></table>

$\gamma_1-\gamma_3$ are used to calculate $h_s$, $\gamma_4-\gamma_6$ are used to calculate $K_s$, $\gamma_7-\gamma_8$ are used to calculate $\lambda$, and $\gamma_9-\gamma_{10}$ are used to calculate $\theta_s$).
<table><caption>**Table S7.** mHM-RE 3 model's global calibration parameters ($\gamma$) related to soil moisture (SM) infiltration process, based on the PTF from Rawls and Brakensiek (1989). Other mHM processes are parameterized similarly to the original mHM model as defined in Table S3.</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td rowspan="25">Infiltration</td><td>γ<sub>1</sub></td><td>PTF-h<sub>S1</sub></td><td>2.67</td><td>10.68</td><td>5.340</td></tr><tr><td>γ<sub>2</sub></td><td>PTF-h<sub>S2</sub></td><td>0.093</td><td>0.370</td><td>0.185</td></tr><tr><td>γ<sub>3</sub></td><td>PTF-h<sub>S3</sub></td><td>1.242</td><td>4.968</td><td>2.484</td></tr><tr><td>γ<sub>4</sub></td><td>PTF-h<sub>S4</sub></td><td>0.001</td><td>0.004</td><td>0.002</td></tr><tr><td>γ<sub>5</sub></td><td>PTF-h<sub>S5</sub></td><td>0.022</td><td>0.088</td><td>0.044</td></tr><tr><td>γ<sub>6</sub></td><td>PTF-h<sub>S6</sub></td><td>0.3085</td><td>1.234</td><td>0.617</td></tr><tr><td>γ<sub>7</sub></td><td>PTF-h<sub>S7</sub></td><td>0.0005</td><td>0.002</td><td>0.0014</td></tr><tr><td>γ<sub>8</sub></td><td>PTF-h<sub>S8</sub></td><td>0.0045</td><td>0.018</td><td>0.009</td></tr><tr><td>γ<sub>9</sub></td><td>PTF-h<sub>S9</sub></td><td>5.0E-6</td><td>20E-6</td><td>12.8E-6</td></tr><tr><td>γ<sub>10</sub></td><td>PTF-h<sub>S10</sub></td><td>0.0045</td><td>0.018</td><td>0.009</td></tr><tr><td>γ<sub>11</sub></td><td>PTF-h<sub>S11</sub></td><td>3.5E-4</td><td>14.0E-4</td><td>7.3E-4</td></tr><tr><td>γ<sub>12</sub></td><td>PTF-h<sub>S12</sub></td><td>2.5E-6</td><td>10.0E-6</td><td>5.4E-6</td></tr><tr><td>γ<sub>13</sub></td><td>PTF-h<sub>S13</sub></td><td>0.25</td><td>1.00</td><td>0.50</td></tr><tr><td>γ<sub>14</sub></td><td>PTF-K<sub>S1</sub></td><td>9.762</td><td>39.047</td><td>19.524</td></tr><tr><td>γ<sub>15</sub></td><td>PTF-K<sub>S2</sub></td><td>-17.937</td><td>-4.484</td><td>-8.969</td></tr><tr><td>γ<sub>16</sub></td><td>PTF-K<sub>S3</sub></td><td>-0.056</td><td>0.014</td><td>-0.028</td></tr><tr><td>γ<sub>17</sub></td><td>PTF-K<sub>S4</sub></td><td>9.1E-5</td><td>3.6E-4</td><td>1.8E-4</td></tr><tr><td>γ<sub>18</sub></td><td>PTF-K<sub>S5</sub></td><td>-0.019</td><td>-0.005</td><td>-0.009</td></tr><tr><td>γ<sub>19</sub></td><td>PTF-K<sub>S6</sub></td><td>-16.790</td><td>-4.198</td><td>-8.395</td></tr><tr><td>γ<sub>20</sub></td><td>PTF-K<sub>S7</sub></td><td>0.039</td><td>0.155</td><td>0.078</td></tr><tr><td>γ<sub>21</sub></td><td>PTF-K<sub>S8</sub></td><td>-5.96E-3</td><td>-1.49E-3</td><td>-2.98E-3</td></tr><tr><td>γ<sub>22</sub></td><td>PTF-K<sub>S9</sub></td><td>-0.039</td><td>-0.010</td><td>-0.020</td></tr><tr><td>γ<sub>23</sub></td><td>PTF-K<sub>S10</sub></td><td>8.65E-6</td><td>34.6E-6</td><td>17.3E-6</td></tr><tr><td>γ<sub>24</sub></td><td>PTF-K<sub>S11</sub></td><td>0.014</td><td>0.055</td><td>0.027</td></tr><tr><td>γ<sub>25</sub></td><td>PTF-K<sub>S12</sub></td><td>0.72E-3</td><td>2.87E-3</td><td>1.43E-3</td></tr><tr><td>γ<sub>26</sub></td><td>PTF-K<sub>S13</sub></td><td>-7.0E-6</td><td>-1.75E-6</td><td>-3.5E-6</td></tr><tr><td>γ<sub>27</sub></td><td>PTF-λ<sub>1</sub></td><td>0.392</td><td>1.568</td><td>0.784</td></tr><tr><td>γ<sub>28</sub></td><td>PTF-λ<sub>2</sub></td><td>0.009</td><td>0.036</td><td>0.018</td></tr><tr><td>γ<sub>29</sub></td><td>PTF-λ<sub>3</sub></td><td>0.531</td><td>2.124</td><td>1.062</td></tr><tr><td>γ<sub>30</sub></td><td>PTF-λ<sub>4</sub></td><td>2.5E-5</td><td>10E-5</td><td>5.3E-5</td></tr><tr><td>γ<sub>31</sub></td><td>PTF-λ<sub>5</sub></td><td>0.002</td><td>0.006</td><td>0.003</td></tr><tr><td>γ<sub>32</sub></td><td>PTF-λ<sub>6</sub></td><td>0.555</td><td>2.222</td><td>1.111</td></tr><tr><td>γ<sub>33</sub></td><td>PTF-λ<sub>7</sub></td><td>0.016</td><td>0.062</td><td>0.031</td></tr><tr><td>γ<sub>34</sub></td><td>PTF-λ<sub>8</sub></td><td>1.5E-4</td><td>6.0E-4</td><td>2.7E-4</td></tr><tr><td>γ<sub>35</sub></td><td>PTF-λ<sub>9</sub></td><td>0.003</td><td>0.012</td><td>0.006</td></tr><tr><td>γ<sub>36</sub></td><td>PTF-λ<sub>10</sub></td><td>1.0E-6</td><td>4.0E-6</td><td>2.0E-6</td></tr><tr><td>γ<sub>37</sub></td><td>PTF-λ<sub>11</sub></td><td>0.004</td><td>0.016</td><td>0.008</td></tr><tr><td>γ<sub>38</sub></td><td>PTF-λ<sub>12</sub></td><td>0.0035</td><td>0.014</td><td>0.007</td></tr><tr><td>γ<sub>39</sub></td><td>PTF-θ<sub>s1</sub></td><td>0.700</td><td>0.999</td><td>0.930</td></tr><tr><td>γ<sub>40</sub></td><td>PTF-θ<sub>s2</sub></td><td>1.00</td><td>5.00</td><td>2.65</td></tr></tbody></table>

June 26, 2025, 11:27pm
<table><caption>**Table S7** continued from previous page</caption><thead><tr><th>Process description</th><th>γ</th><th>Parameter name</th><th>Min.</th><th>Max.</th><th>Init.</th></tr></thead><tbody><tr><td>γ<sub>41</sub></td><td>PTF-θ<sub>r1</sub></td><td></td><td>0.009</td><td>0.036</td><td>0.018</td></tr><tr><td>γ<sub>42</sub></td><td>PTF-θ<sub>r2</sub></td><td></td><td>4.5E-4</td><td>18.0E-4</td><td>8.7E-4</td></tr><tr><td>γ<sub>43</sub></td><td>PTF-θ<sub>r3</sub></td><td></td><td>0.0025</td><td>0.010</td><td>0.005</td></tr><tr><td>γ<sub>44</sub></td><td>PTF-θ<sub>r4</sub></td><td></td><td>0.015</td><td>0.058</td><td>0.029</td></tr><tr><td>γ<sub>45</sub></td><td>PTF-θ<sub>r5</sub></td><td></td><td>1.0E-4</td><td>4.0E-4</td><td>2.0E-4</td></tr><tr><td>γ<sub>46</sub></td><td>PTF-θ<sub>r6</sub></td><td></td><td>0.0005</td><td>0.002</td><td>0.001</td></tr><tr><td>γ<sub>47</sub></td><td>PTF-θ<sub>r7</sub></td><td></td><td>1.0E-4</td><td>4.0E-4</td><td>2.0E-4</td></tr><tr><td>γ<sub>48</sub></td><td>PTF-θ<sub>r8</sub></td><td></td><td>1.5E-4</td><td>6.0E-4</td><td>3.1E-4</td></tr><tr><td>γ<sub>49</sub></td><td>PTF-θ<sub>r9</sub></td><td></td><td>0.001</td><td>0.004</td><td>0.002</td></tr></tbody></table>

γ₁-γ₁₃ are used to calculate *h*<sub>s</sub>, γ₁₄-γ₂₆ are used to calculate *K*<sub>s</sub>, γ₂₇-γ₃₈ are used to calculate λ, γ₃₉-γ₄₀ are used to calculate θ<sub>s</sub>, and γ₄₁-γ₄₉ are used to calculate θ<sub>r</sub>).
<table><caption>**Table S8.** Soil hydraulic parameters for the 11 USDA soil classes estimated from Cosby et al. (1984) PTF (a index) and **mHM-RE 1** model (b index).</caption><thead><tr><th>Texture class</th><th>Sand</th><th>Clay</th><th>Silt</th><th>λ<sup>a</sup></th><th>λ<sup>b</sup></th><th>h<sub>s</sub><sup>a</sup></th><th>h<sub>s</sub><sup>b</sup></th><th>θ<sub>s</sub><sup>a</sup></th><th>θ<sub>s</sub><sup>b</sup></th><th>K<sub>s</sub><sup>a</sup></th><th>K<sub>s</sub><sup>b</sup></th></tr></thead><tbody><tr><td>Sand</td><td>92</td><td>3</td><td>5</td><td>0.36</td><td>N.A.</td><td>6.92</td><td>N.A.</td><td>0.34</td><td>N.A.</td><td>16.78</td><td>N.A.</td></tr><tr><td>Loamy sand</td><td>82</td><td>6</td><td>12</td><td>0.23</td><td>0.24</td><td>3.63</td><td>12.59</td><td>0.42</td><td>0.41</td><td>5.07</td><td>3.08</td></tr><tr><td>Sandy loam</td><td>58</td><td>10</td><td>32</td><td>0.21</td><td>0.20</td><td>14.13</td><td>15.85</td><td>0.43</td><td>0.37</td><td>1.88</td><td>1.79</td></tr><tr><td>Loam</td><td>43</td><td>18</td><td>39</td><td>0.19</td><td>0.16</td><td>35.48</td><td>41.69</td><td>0.44</td><td>0.40</td><td>1.22</td><td>1.06</td></tr><tr><td>Silt loam</td><td>17</td><td>13</td><td>70</td><td>0.19</td><td>0.17</td><td>75.86</td><td>48.05</td><td>0.48</td><td>0.50</td><td>1.01</td><td>0.69</td></tr><tr><td>Sandy clay loam</td><td>58</td><td>27</td><td>15</td><td>0.15</td><td>0.23</td><td>13.49</td><td>10.00</td><td>0.40</td><td>0.39</td><td>1.60</td><td>1.48</td></tr><tr><td>Clay loam</td><td>32</td><td>34</td><td>34</td><td>0.12</td><td>0.13</td><td>26.30</td><td>44.69</td><td>0.47</td><td>0.42</td><td>0.88</td><td>0.72</td></tr><tr><td>Silty clay loam</td><td>10</td><td>34</td><td>56</td><td>0.11</td><td>0.11</td><td>61.66</td><td>49.35</td><td>0.46</td><td>0.50</td><td>0.73</td><td>0.48</td></tr><tr><td>Sandy clay</td><td>52</td><td>42</td><td>6</td><td>0.09</td><td>N.A.</td><td>9.77</td><td>N.A.</td><td>0.41</td><td>N.A.</td><td>2.60</td><td>N.A.</td></tr><tr><td>Silty clay</td><td>6</td><td>47</td><td>47</td><td>0.10</td><td>N.A.</td><td>32.36</td><td>N.A.</td><td>0.47</td><td>N.A.</td><td>0.48</td><td>N.A.</td></tr><tr><td>Clay</td><td>22</td><td>58</td><td>20</td><td>0.09</td><td>N.A.</td><td>46.77</td><td>N.A.</td><td>0.47</td><td>N.A.</td><td>0.35</td><td>N.A.</td></tr></tbody></table>

Note that the soil dataset for the 201 basins analyzed in this study does not include the USDA textural classes of Sand, Sandy clay, Silty clay, and Slay during the calibration. λ unit is [-], h_s unit is [cm], θ_s unit is [cm³ cm⁻³], K_s unit is [cm hour⁻¹], and Sand, Clay, and Silt are in percentage.

<table><caption>**Table S9.** Similar to Table S8, but the data were estimated from Campbell (1985) PTF (a index) and **mHM-RE 2** model (b index).</caption><thead><tr><th>Texture class</th><th>Sand</th><th>Clay</th><th>Silt</th><th>λ<sup>a</sup></th><th>λ<sup>b</sup></th><th>h<sub>s</sub><sup>a</sup></th><th>h<sub>s</sub><sup>b</sup></th><th>θ<sub>s</sub><sup>a</sup></th><th>θ<sub>s</sub><sup>b</sup></th><th>K<sub>s</sub><sup>a</sup></th><th>K<sub>s</sub><sup>b</sup></th></tr></thead><tbody><tr><td>Sand</td><td>92</td><td>3</td><td>5</td><td>0.48</td><td>N.A.</td><td>9.49</td><td>N.A.</td><td>0.42</td><td>N.A.</td><td>3.27</td><td>N.A.</td></tr><tr><td>Loamy sand</td><td>82</td><td>6</td><td>12</td><td>0.36</td><td>0.30</td><td>11.72</td><td>14.91</td><td>0.43</td><td>0.41</td><td>1.21</td><td>0.63</td></tr><tr><td>Sandy loam</td><td>58</td><td>10</td><td>32</td><td>0.23</td><td>0.26</td><td>18.48</td><td>11.74</td><td>0.44</td><td>0.49</td><td>0.15</td><td>0.018</td></tr><tr><td>Loam</td><td>43</td><td>18</td><td>39</td><td>0.16</td><td>0.20</td><td>27.75</td><td>29.73</td><td>0.45</td><td>0.49</td><td>0.02</td><td>0.0065</td></tr><tr><td>Silt loam</td><td>17</td><td>13</td><td>70</td><td>0.13</td><td>0.15</td><td>37.48</td><td>28.11</td><td>0.53</td><td>0.51</td><td>0.00651</td><td>0.0057</td></tr><tr><td>Sandy clay loam</td><td>58</td><td>27</td><td>15</td><td>0.15</td><td>0.23</td><td>26.61</td><td>15.66</td><td>0.47</td><td>0.46</td><td>0.02952</td><td>0.014</td></tr><tr><td>Clay loam</td><td>32</td><td>34</td><td>34</td><td>0.10</td><td>0.13</td><td>46.75</td><td>47.60</td><td>0.48</td><td>0.45</td><td>0.00243</td><td>0.0073</td></tr><tr><td>Silty clay loam</td><td>10</td><td>34</td><td>56</td><td>0.08</td><td>0.11</td><td>66.20</td><td>59.13</td><td>0.48</td><td>0.49</td><td>0.00052</td><td>0.00100</td></tr><tr><td>Sandy clay</td><td>52</td><td>42</td><td>6</td><td>0.11</td><td>N.A.</td><td>40.58</td><td>N.A.</td><td>0.50</td><td>N.A.</td><td>0.00453</td><td>N.A.</td></tr><tr><td>Silty clay</td><td>6</td><td>47</td><td>47</td><td>0.05</td><td>N.A.</td><td>93.92</td><td>N.A.</td><td>0.51</td><td>N.A.</td><td>0.00011</td><td>N.A.</td></tr><tr><td>Clay</td><td>22</td><td>58</td><td>20</td><td>0.05</td><td>N.A.</td><td>92.89</td><td>N.A.</td><td>0.51</td><td>N.A.</td><td>0.00012</td><td>N.A.</td></tr></tbody></table>

Note that the soil dataset for the 201 basins analyzed in this study does not include the USDA textural classes of Sand, Sandy clay, Silty clay, and Slay during the calibration. λ unit is [-], h_s unit is [cm], θ_s unit is [cm³ cm⁻³], K_s unit is [cm hour⁻¹], and Sand, Clay, and Silt are in percentage.
<table><caption>**Table S10.** Similar to Table S8, but the data were estimated from Rawls and Brakensiek (1989) PTF (a index) and mHM-RE 3 model (b index).</caption><thead><tr><th>Texture class</th><th>Sand</th><th>Clay</th><th>Silt</th><th>λ<sup>a</sup></th><th>λ<sup>b</sup></th><th>h<sub>s</sub><sup>a</sup></th><th>h<sub>s</sub><sup>b</sup></th><th>θ<sub>s</sub><sup>a</sup></th><th>θ<sub>s</sub><sup>b</sup></th><th>K<sub>s</sub><sup>a</sup></th><th>K<sub>s</sub><sup>b</sup></th><th>θ<sub>r</sub><sup>a</sup></th><th>θ<sub>r</sub><sup>b</sup></th></tr></thead><tbody><tr><td>Sand</td><td>92</td><td>3</td><td>5</td><td>0.69</td><td>N.A.</td><td>15.98</td><td>N.A.</td><td>0.44</td><td>N.A.</td><td>21.00</td><td>N.A.</td><td>0.020</td><td>N.A.</td></tr><tr><td>Loamy sand</td><td>82</td><td>6</td><td>12</td><td>0.55</td><td>0.40</td><td>20.58</td><td>30.34</td><td>0.44</td><td>0.41</td><td>6.11</td><td>13.18</td><td>0.035</td><td>0.033</td></tr><tr><td>Sandy loam</td><td>58</td><td>10</td><td>32</td><td>0.38</td><td>0.33</td><td>30.20</td><td>32.32</td><td>0.45</td><td>0.45</td><td>2.59</td><td>1.58</td><td>0.041</td><td>0.045</td></tr><tr><td>Loam</td><td>43</td><td>18</td><td>39</td><td>0.25</td><td>0.19</td><td>40.12</td><td>37.22</td><td>0.46</td><td>0.46</td><td>0.68</td><td>1.12</td><td>0.027</td><td>0.050</td></tr><tr><td>Silt loam</td><td>17</td><td>13</td><td>70</td><td>0.23</td><td>0.24</td><td>50.87</td><td>39.07</td><td>0.50</td><td>0.49</td><td>1.32</td><td>0.27</td><td>0.015</td><td>0.047</td></tr><tr><td>Sandy clay loam</td><td>58</td><td>27</td><td>15</td><td>0.32</td><td>0.56</td><td>59.41</td><td>49.99</td><td>0.40</td><td>0.48</td><td>0.43</td><td>0.29</td><td>0.068</td><td>0.025</td></tr><tr><td>Clay loam</td><td>32</td><td>34</td><td>34</td><td>0.24</td><td>0.18</td><td>56.43</td><td>42.67</td><td>0.46</td><td>0.41</td><td>0.23</td><td>0.20</td><td>0.075</td><td>0.042</td></tr><tr><td>Silty clay loam</td><td>10</td><td>34</td><td>56</td><td>0.18</td><td>0.14</td><td>70.33</td><td>47.27</td><td>0.47</td><td>0.42</td><td>0.15</td><td>0.12</td><td>0.040</td><td>0.053</td></tr><tr><td>Sandy clay</td><td>52</td><td>42</td><td>6</td><td>0.22</td><td>N.A.</td><td>79.48</td><td>N.A.</td><td>0.43</td><td>N.A.</td><td>0.12</td><td>N.A.</td><td>0.109</td><td>N.A.</td></tr><tr><td>Silty clay</td><td>6</td><td>47</td><td>47</td><td>0.15</td><td>N.A.</td><td>76.54</td><td>N.A.</td><td>0.48</td><td>N.A.</td><td>0.09</td><td>N.A.</td><td>0.056</td><td>N.A.</td></tr><tr><td>Clay</td><td>22</td><td>58</td><td>20</td><td>0.17</td><td>N.A.</td><td>85.60</td><td>N.A.</td><td>0.48</td><td>N.A.</td><td>0.06</td><td>N.A.</td><td>0.090</td><td>N.A.</td></tr></tbody></table>

Note that the soil dataset for the 201 basins analyzed in this study does not include the USDA textural classes of Sand, Sandy clay, Silty clay, and Slay during the calibration. λ unit is [-], h_s unit is [cm], θ_s unit is [cm³ cm⁻³], K_s unit is [cm hour⁻¹], and Sand, Clay, and Silt are in percentage.
References

Boeing, F., Rakovec, O., Kumar, R., Samaniego, L., Schrön, M., Hildebrandt, A., ... Marx, A. (2022). High-resolution drought simulations and comparison to soil moisture observations in Germany. *Hydrology and Earth System Sciences*, 26(19), 5137–5161. doi: 10.5194/hess-26-5137-2022

Campbell, G. S. (1985). *Soil physics with BASIC: transport models for soil-plant systems* (No. 14). Amsterdam New York: Elsevier.

Cosby, B. J., Hornberger, G. M., Clapp, R. B., & Ginn, T. R. (1984). A Statistical Exploration of the Relationships of Soil Moisture Characteristics to the Physical Properties of Soils. *Water Resources Research*, 20(6), 682–690. doi: 10.1029/WR020i006p00682

EEA. (2020). *CORINE Land Cover 1990, 2000, 2006, 2012, 2018 vector/raster 100 m.* Copernicus Land Monitoring Service. Retrieved 2024-10-21, from https://land.copernicus.eu/en/products/corine-land-cover

Hargreaves, G. H., & Samani, Z. A. (1985). Reference Crop Evapotranspiration from Temperature. *Applied Engineering in Agriculture*, 1(2), 96–99. doi: 10.13031/2013.26773

Hartmann, J., & Moosdorf, N. (2012a). *Global Lithological Map Database v1.0 (gridded to 0.5° spatial resolution)*. PANGAEA. doi: https://doi.org/10.1594/PANGAEA.788537

Hartmann, J., & Moosdorf, N. (2012b). The new global lithological map database GLiM: A representation of rock properties at the Earth surface. *Geochemistry, Geophysics, Geosystems*, 13(12). doi: 10.1029/2012GC004370

Jarvis, A., Reuter, H. I., Nelson, A., & Guevara, E. (2008). *Hole-filled seamless SRTM data V4*. International Centre for Tropical Agriculture (CIAT). Retrieved 2025-06-26, from https://srtm.csi.cgiar.org
Myneni, R., Knyazikhin, Yuri, & Park, T. (2015). *Myd15a2h modis/aqua leaf area index/fpar 8-day l4 global 500m sin grid. nasa lp daac*. Boston University and MODAPS SIPS - NASA. doi: http://doi.org/10.5067/MODIS/MYD15A2H.006

Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Y.-W., ... Papale, D. (2020). The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. *Scientific Data*, 7(1), 225. doi: 10.1038/s41597-020-0534-3

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., & Rossiter, D. (2021). SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. *SOIL*, 7(1), 217-240. doi: 10.5194/soil-7-217-2021

Rawls, W. J., & Brakensiek, D. L. (1989). Estimation of Soil Water Retention and Hydraulic Properties. In H. J. Morel-Seytoux (Ed.), *Unsaturated Flow in Hydrologic Modeling* (pp. 275-300). Dordrecht: Springer Netherlands. doi: 10.1007/978-94-009-2352-2_10

Samaniego, L., Kumar, R., & Jackisch, C. (2011). Predictions in a data-sparse region using a regionalized grid-based hydrologic model driven by remotely sensed data. *Hydrology Research*, 42(5), 338-355. doi: 10.2166/nh.2011.156

Zacharias, S., Loescher, H. W., Bogena, H., Kiese, R., Schrön, M., Attinger, S., ... Vereecken, H. (2024). Fifteen Years of Integrated Terrestrial Environmental Observatories (TERENO) in Germany: Functions, Services, and Lessons Learned. *Earth's Future*, 12(6). doi: 10.1029/2024EF004510