MPR 1.0: A stand-alone Multiscale Parameter Regionalization Tool for Improved Parameter Estimation of Land Surface Models

Distributed environmental models such as land surface models (LSM) require model parameters in each spatial modelling unit (e.g. grid cell), thereby leading to a high-dimensional parameter space. One approach to decrease the dimensionality of parameter space in these models is to use regularization techniques. One such highly efficient technique is the Multiscale Parameter Regionalization (MPR) framework that translates high-resolution predictor variables (e.g., soil textural properties) into model parameters (e.g., porosity) via transfer functions (TFs) and upscaling operators that are suitable for every 5 modeled process. This framework yields seamless model parameters at multiple scales and locations in an effective manner. However, integration of MPR into existing modeling workflows has been hindered thus far by hard-coded configurations and non-modular software designs. For these reasons, we redesigned MPR as a model-agnostic, stand-alone tool. It is a useful software for creating graphs of netCDF variables, wherein each node is a variable and the links consist of TFs and/or upscaling operators. In this study, we present and verify our tool against a previous version, which was implemented in the mesoscale 10 hydrologic model mHM (www.ufz.de/mhm). By using this tool for the generation of continental-scale soil hydraulic parameters applicable to different models (Noah-MP and HTESSEL), we showcase its general functionality and flexibility. Further, using model parameters estimated by the MPR tool leads to significant changes in long-term estimates of evapotranspiration, as compared to their default parameterizations. For example, a change of up to 25% in long-term evapotranspiration flux is observed in Noah-MP and HTESSEL in the Mississippi River basin. We postulate that use of the stand-alone MPR tool will 15 considerably increase the transparency and reproducibility of the parameter estimation process in distributed (environmental) models. It will also allow a rigorous uncertainty estimation related to the errors of the predictors (e.g., soil texture fields), transfer function and its parameters, and remapping (or upscaling) algorithms. 1 https://doi.org/10.5194/gmd-2021-103 Preprint. Discussion started: 26 May 2021 c © Author(s) 2021. CC BY 4.0 License.

1 Introduction 20 Distributed environmental models simulate key fluxes and states of the atmosphere, land surface, and subsurface for a given spatial domain and time period (e.g., CLM (Andre et al., 2020), JULES (Best et al., 2011), or ORCIDEE (Krinner et al., 2005). The underlying physical processes are simplified with parameterizations that are manifested as computer algorithms.
Parameterizations are idealized representations of reality and as such there is inherent uncertainty in their formulation. They require additional variables and model parameters in order to perform simulations. The latter could be constant, or it could 25 be spatially and temporally variable over the simulation domain (i.e., the so-called distributed model parameters). Constant model parameters do not allow accurate characterization of environmental processes over a range of climatic regimes and geo-physical properties (Samaniego et al., 2010;Beck et al., 2016). The number of distributed model parameters tend to scale linearly with the number of spatio-temporal units, which is defined by the coordinates along each dimension of the parameterized process. Model parameters are often fine-tuned to match observed fluxes and states of physical processes in 30 various fields, such as hydrology (Duan et al., 1992;Gupta et al., 1998;Zink et al., 2017;Pagliero et al., 2019), Earth System Sciences (Sellers et al., 1989;Troy et al., 2008), and hydraulics (Shoarinezhad et al., 2020). Currently, there exists a plethora of methods for parameter estimation , which is the process of estimating a set of model parameters over the whole domain and their respective distribution functions. Such methods include classification linked through lookup tables (for example, ECMWF, 2019; Andre et al., 2020), direct calibration (for example, Li et al., 2018;Arheimer et al., 35 Providing a library of remapping schemes is crucial, as the parameterizations in environmental models are not applied on the scales at which they were derived. For example, the Richardson & Richards' equation (Richardson, 1922;Richards, 90 1931) describing unsaturated water flow through porous media at the representative elemental volume scale is often used at the landscape scale (say 10 2 km). The inherent uncertainty of the physical parameters describing this phenomenon needs to be adequately considered by the choice of transfer functions and upscaling operators (Montzka et al., 2017;Vereecken et al., 2019). More generally, flow rates such as saturated hydraulic conductivity are common parameters in environmental models, and their scaling behavior should be considered (Zhu and Mohanty, 2002;Kumar et al., 2013b). Additionally, their anisotropic 95 properties necessitate a dimension-dependent selection of upscaling operators (e.g., harmonic and arithmetic mean).
In this paper, we first present a working example to illustrate current challenges in the estimation of distributed parameter fields in environmental models. We then outline different ways to address this task and demonstrate how the new generic and agnostic MPR tool is designed and configured to meet the requirements. The section elucidating the configuration of MPR is followed by a detailed description of how to interface MPR through a stand-alone executable, as well as through its API 100 (application programming interface). We demonstrate its tight coupling to the hydrologic model mHM and its capabilities with regard to the reproduction of the original mHM model behavior. Furthermore, we demonstrate the versatility of the MPR tool by reproducing an open source dataset (EU-SoilHydroGrids (Tóth et al., 2017)) containing soil hydraulic properties derived from a set of TFs. The effects on long-term evapotranspiration from the coupling of MPR to state-of-the-art land surface models (Noah-MP and HTESSEL) are also shown for an effective TF application and developing a regridding tool applicable to any 105 model that requires distributed parameters (e.g., LSM or environmental models) 2 A minimal working example in environmental modelling

Objective
For demonstration purposes we define an objective that is commonly encountered in environmental modeling. For a hypothetical inter-comparison project, the influence of different parameter estimation schemes for soil parameters is investigated for a 110 given domain along with two different resolutions and three different model-specific grid types.
A common parameter present in many environmental models is soil porosity (θ s ), which denotes the pore volume fraction of the total soil volume in the vadose zone. The SoilGrids (SG) dataset  provides soil physical properties at a high resolution ( 1 480 • ). From the extensive literature on pedotransfer functions (TFs for soil parameters) (Patil and Singh, 2016), we selected a TF for estimating θ s based on bulk density, organic matter, clay, and sand content (Weynants et al.,115 2009). The southeastern United States, which includes the state of Florida, was chosen as the domain of interest due to high heterogeneity in the physical properties of the soil in this region. We selected different grid layouts that are often used in different modeling disciplines. These are regular rectangular grids generally used in distributed environmental modeling (Ma et al., 2017;Zink et al., 2017), icosahedral grids representing the group of geodesic grids increasingly used in the Earth System Science community (Zängl et al., 2015;Skamarock et al., 2012), and polygons or hydrologic response units (HRUs) often used 120 in hydrology or the soil sciences (Wellen et al., 2015). The exemplary workflow depicted in Fig. 1 shows the high-resolution predictor variables (shown in the lower left corner) that are passed to the transfer function (TF) to derive the resulting field porosity (lower right corner) at the predictor resolution.
It exhibits considerable heterogeneity in various gradients between the southeast and northwest regions of the domain, coastal and inland areas, riverbeds, and mountainous areas. The target variable is then upscaled to six different spatial grids at different The lower grid is the same as the one used for the North American Land Data Assimilation System 2 forcing dataset (Mitchell et al., 2004). The latter dataset covers the domain of the conterminous United States (CONUS) and has been used for different LSMs (Xia et al., 2012). The two icosahedral grids specified by identifier R02B04 and R02B05 ( Fig. 1 the center column in the top panel (Zängl et al., 2015)) can be used for the configuration of LSM JSBACH (Mauritsen et al., 2019). Finally, for HRU-based (Flügel, 1995) models such as SWAT (P. W. Gassman et al., 2007) or PRMS (Leavesley et al., 1983). The This parameter estimation routine, which is followed by simulations of the distributed model, might then be used iteratively by a calibration routine in order to optimize TF parameters. These are global parameters of the TF whose application to the predictors is used to derive the effective distributed model parameters. Next, a workflow to derive these parameters in an 140 efficient and consistent manner is presented.

Options for parameter estimation workflow
There exists a range of different options for the workflow described in the previous section, which we briefly describe here.
A few existing software tools and GCM couplers can be used to perform the two key steps of applying a TF and remapping a multi-dimensional grid onto an another: monotonically ordered integers and can be shared among multiple variables. Each variable and dimension has a name, and the referencing of dimensions is done using these names. and Upscalers are designed as arrays of respective objects whose specific properties are set using the correct list index.
The term Data_Arrays refers to any n-dimensional variable and serves as a generic term for predictors and target variables. In lines 1ff. in Fig. 2, there are four Data_Arrays specified: the first three are read directly from the file, while the 215 fourth is calculated from the former. The property from_data_arrays specifies the array of predictors to be used. They also appear in the TF equation that is supplied as a string in transfer_func. It follows the Fortran syntax for operators, brackets, and some elemental mathematical functions (see Appendix Table E1 for a list of possible operators). Users can use any parameter defined in Parameters in a TF. The TF string is not dynamically evaluated during execution because it leads to unnecessarily high computational run times. Additionally, the effort required to modify the Fortran source code each time a 220 new function is used would substantially diminish the user-friendliness of the MPR tool. Instead, a Python preprocessor script is implemented that interprets and adds the transfer function from the namelist and modifies the source code accordingly. At runtime, simple search and replacement routines translate the string into a unique function ID that is checked against all unique function IDs in the source code.
The target coordinates and associated upscaling operators are set with target_coord_names and upscale_ops.
seperate file containing only the parameters subject to optimization. This file is optional and its parameters replace previously read parameters in the case of duplicates.

235
The target coordinates are specified in the Coordinates section (l. 26ff. in Fig. 2). The user explicitly specifies the boundaries of the soil layers by values in this example (coord_from_values). These refer to the stagger of each cell or horizon (coord_stagger). In this case, the first cell does not have a neighboring cell boundary, and as such the user must specify the coordinate boundary to provide a start value (coord_from_values_bound). A two-dimensional coordinate serving as the target grid is read from the file (coord_from_file). Its associated dimensions are specified through coord_sub_dims.

240
Alternatively, coordinate values can be specified from a range if the user provides values for start, step, and count.
The Main section (l. 37 ff. in Fig. 2) configures general information. The link from the target coordinates to the respective source coordinates during upscaling is constructed through the coordinate_groups entries. During upscaling, each source coordinate needs to receive a target coordinate. If the source coordinate shares the same coordinate_group as one of the target coordinates, upscaling is performed for that coordinate pair. Users can enter an arbitrary number of groups, although a 245 dimension system with X, Y, Z, and T is most commonly defined. Finally, the path to the output file is set.
A more exhaustive description of the aforementioned configuration can be found in online documentation.

Interdependence of parameters
Another example configuration highlighting the capability of MPR for creation of netCDF variable graphs is shown in Fig. 3,where each node is a Data_Array and the links consisting of TFs and/or upscaling operations. It visualizes the dependency graph 250 for two of the parameters required for the mHM in its standard configuration. The blue ellipses denote the predictor variables used. While land_cover is a 3-dimensional array with coordinates, year, latitude, longitude (t1, y, x), lai_class is a 3-dimensional array with coordinates, month of year, latitude, and longitude (t2, y, x). The TF for the model parameter PET_LAI_corr_factor requires both LAI and land cover information. Both predictors need to be broadcast to a 4-dimensional array (t1, t2, y, x), and so temporary arrays are created. The TF for the model parameter Aerodyn_resist 255 requires information on canopy and wind measurement height, while the wind measurement height is derived from the former alone. For the calculation of canopy_height, a max-normalized LAI is needed for each cell. The Data_Array contains this information, and then its broadcast variant LAI_max_t2_4D will be the same shape as the other predictors of the TF responsible for producing canopy_height. The two model parameters highlighted by the red ellipsis are finally upscaled to the target model resolution with coordinates modeling year, month of year, low resolution latitude, and low resolution longitude.

Interfacing MPR library
While the previous section introduced the use of MPR as a standalone wrapper, we anticipate a tight coupling of the MPR code to the main modeling code. We intend to impose maximum reusability of the MPR API and ease its implementation into other libraries. Top-level objects such as Data_Arrays or Coordinates can be reused multiple times and can also be written to and read from the disk as requested. The Fortran API is used here for this purpose, which is based on the object-oriented programming 265 paradigm and exposes four main objects (derived type in Fortran), which are also present in the namelist configuration (see  previous section): Coordinate, Parameter, Data_Array and CoordUpscaler. After initialization, these types are stored in global arrays which allow for cross-referencing with other types. They are briefly described in this section and additional detail can be found in the online documentation.

Coordinate type 270
One central derived type in MPR is the Coordinate. First and foremost, it stores the boundaries of its n cells. In the case of a (one-dimensional) coordinate variable, each cell has two boundaries (v = 2) and as cells are contiguous, the boundaries do not need to be stored in a (v, n)-shaped array but in an (n + 1) array. This property is termed as boundaries1d. In the case of a two-dimensional auxiliary coordinate variable, the number of boundaries/edges can vary. To obtain a general formulation, we stored the boundaries of both dimensions in an (v, n)-shaped array. These properties are termed boundaries2dDim1 275 and boundaries2dDim2. Additionally, the cell centers are stored in centers2dDim1 and centers2dDim2.

Parameter type
Parameters are objects with names and numerical values assigned to them.

Data_Array type
The main derived type is the Data_Array, which can be read directly from an existing netCDF file or computed from other 280 Data_Arrays and/or Parameters using TF and upscaling operators. It stores multidimensional data which must be passed to multiple other routines. In order to have a common sparse data container, non-masked cells are stored in a one-dimensional array data with type real, regardless of the number of underlying dimensions. A Boolean mask of the data is stored in a flattened, one-dimensional array, reshapedMask. A one-dimensional array pointing to its associated Coordinates is set as coords. This holds the shape information of the original uncompressed data. These core properties, reshapedMask, 285 data, and coords, are pointed at from within the wrapper type InputFieldContainer. It is used for referencing the core properties of the Data_Array and are usually passed as an argument for the TransferHelper type, which is described in the next section.

Application of Transfer Functions
The type TransferHelper is intended to be used for the initialization of Data_Array, which originates from either a file 290 or from other Data_Arrays via TFs. During its initialization, it performs checks on the passed InputFieldContainer, checks their Coordinates, optionally concatenates Data_Arrays, and optionally adapts their masks. After these checks, the TF is applied. TFs are designed as a subroutine accepting an arbitrarily-sized list of InputFieldContainer (its data property is accessed only) and Parameters, respectively. An abstract interface for TFs thus allows for a variable number of predictor arrays and parameters. A one-dimensional array is returned. We provide a template for a subroutine containing a TF 295 so that users can modify the source code and implement their own TFs (see Appendix Fig. D1). This enables the integration of more complex mathematical formulations such as artificial neural networks or support vector machines. Based on the majority of TFs used in existing LSMs (Van Looy et al., 2017), we anticipate users primarily using TFs that employ common mathematical operators such as multiplication and division. Such equations can be automatically parsed from configuration files and inserted into the Fortran source code using a Python preprocessor, which is described next.

300
Each TF string occurring in a namelist is translated into a unique identifier. For example, string exp(a5 + c5 * sand + d5 * bd + e5 * om) / unit_conversion (Fig. 2) can be set in the configuration file and will then be analyzed and processed by a parser routine. MPR replaces all parameters (a5, c5, d5, e5, unit_conversion), and variable names (sand, bd, om) by identifiers, translating the resulting string (exp (p1 + p2 * x1 + p3 * x2 + p4 * x3) / p5) into a unique function identifier (exp_bs_p1_pl_p2_ti_x1_pl_p3_ti_x2_pl_p4_ti_x3_be_di_p5). This identifier represents the exact mathematical function that can be used for multiple applications with different Data_Arrays and Parameters. The number of TFs contained in the source code is thereby reduced, and duplications of TFs are avoided.

API for Upscaling
For the upscaling step, two Coordinates sharing the same dimensions are compared. For each cell of the target grid, the underlying n source cells (subcells) are stored in an ids array and their relative contributing area is stored in a weights array. These properties are contained in the type CoordUpscaler. They can be initialized from existing remapping weights stored in the netCDF file format following the SCRIP convention (Jones, 2010). By default, first-order conservative remapping 315 is used for weight calculation, so that the integral of the values is preserved in the presence of missing values. Two-dimensional auxiliary coordinate variables are mapped using a simple algorithm that only checks that the cell center of the source cell is within the target cell boundaries and assigns equal weights to each source cell. In the future, more sophisticated remapping schemes for two-dimensional variables will be implemented.
The upscaling of a Data_Array is conducted through the wrapper type UpscaleHelper, which is also a property of the 320 Data_Array type. It consists of a one-dimensional array pointing to the source coordinates sourceCoords and target coordinates targetCoords, as well as the upscale operator names for each target coordinate upscaleOperatorNames. The Upscaler type performs multiple checks on the source and target coordinates of Data_Array. If applicable, it transposes or broadcasts the Data_Array if the order of source and target coordinates does not match. Multiple CoordUpscaler instances can be combined with an aggregated CoordUpscaler object, effectively combining the weights and subcell IDs 325 when the user specifies the same upscaling operator for multiple coordinates of the target grid. The upscaling operation is then executed separately for each group of target dimensions using the same upscaling operator. There are a number of standard upscaling operators, such as the minimum, maximum, or weighted generalized mean function, which employ the power parameter (see Appendix Table E2). Users can easily add another upscaling function to the source code as long as it effectively aggregates the variability of the subcells into one value (see Appendix Fig. D2).

New features of current MPR release
The current set of features encompasses the functionality of the previous version of MPR as used in the mHM source code (Samaniego et al., 2019). This implementation lacks some key features that are highlighted by a comparison of Fig. 1 with Fig. 1 in Kumar et al. (2013a). First, the new MPR library is modularized and refactored, so it does not depend on mHM and its integrated optimization algorithms. Second, the original mHM implementation required grids to be rectangular, and the 335 target resolution was a multiple of the source resolution. Finally, predictors, TFs, upscaling operators, and target grids can now be freely chosen and recombined in an MPR configuration file.
Furthermore, it allows for a modular and flexible configuration of parameter estimation without common preprocessing steps of the input files (resorting coordinates, renaming variables, adding missing meta information, applying unit conversion, etc.).
Use of MPR is easy and intuitive, especially in the formulation of TFs, defining new coordinate variables, and assigning them to 340 variables. The initialization of coordinate variables can be performed dynamically depending on existing coordinate variables (e.g., by using its bounds). The user can freely and flexibly enter simple regression-based TFs in a semi-automatic manner.
More sophisticated functions, such as artificial neural networks or support vector machines, can be easily coupled to the code.
The same holds true for upscaling operators. We allowed MPR to be easily integrated in workflows (e.g., auto-calibration) or have an API called from an external code.

Reproducing EU-SoilHydroGrids dataset with MPR
We selected the EU-SoilHydroGrids (SHG) dataset (Tóth et al., 2017) as it is relevant to the Earth System Modeling community and is publicly available. The same holds true for its predictors and the TFs used. We reproduced the aforementioned dataset and showed the salient seamless spatial scaling feature of the MPR tool. The dataset description and its configuration can be found in the Appendices B2, C2 and the MPR configuration in the Supplements.
365 Fig. 4 shows the spatial distribution of K s (soil saturated hydraulic conductivity) on a logarithmic scale for different resolutions and data sources at a depth of 15cm with values ranging across three orders of magnitude from 10 −6 to 10 −3 m s −1 . Each row of the plot shows grids with the same spatial resolution, and each column shows the data sources and processing schemes.
obtained a bias of up to 10% in comparison to a). An investigation of the project webpage revealed that the SoilGrids (SG) dataset is subject to frequent updates depending on the arrival of new source soil profiles. We could not verify the exact version of SG used to create the SHG, but we assume that it was different than the one we used. The high bias observed is a result of the decision tree, as small deviations in predictor values might result in use of a different branch. An additional source for the observed bias could be the projected coordinate system used for the SHG, in contrast to SG being available in the geographical 375 WGS84 coordinate system.
When using MPR to apply the TFs to the SG dataset (Fig. 4c), we obtained a mean relative bias below 0.06 % compared to b), which was only due to rounding errors of the global parameters. SHG and SG are also available at a 1 km resolution, where variables are aggregated using the Geospatial Data Abstraction Library (contributors, 2019) average method (Fig. 4 d) and e)).

380
In accordance with the MPR framework, we derived the 1 km K s field by averaging the high-resolution SG data at 250 m ( Fig. 4c). Because the selected TF is based on a decision tree algorithm, there exist only some discrete values in Fig. 4d) and e), whereas spatial averaging increases the variability of the derived parameter field. The capability of MPR in flexible aggregation of the data is shown in Fig. 4 (g)). A resolution of 25 km was arbitrarily selected as the target resolution, as specific resolutions (e.g., 1 km) are typically not needed by users. By a mere change of one number in the configuration file, it is possible to 385 scale the variable at interest to every user-defined resolution. This analysis highlights the capability of MPR in reproduction of environmental parameters and generation of outputs that meet specific user requirements.

Application with land surface model Noah-MP
We used the land surface model (LSM) Noah-MP (Niu et al., 2011) as an example to showcase the capability of MPR in coupling with state-of-the-art distributed environmental models. The model description and configuration can be found in the 390 Appendices B3, C3, the MPR configuration in the Supplement. We found that the inherent parameter uncertainty that occurs when choosing a TF and an upscaling operator for the soil parameters of the model eventually leads to differences in long-term mean annual evapotranspiration (ET) flux, up to 20% (Fig. 5) in relation to its default setup. differences (in %) of the field with respect to the default simulation in row 1. The relative differences observed when using the MPR approach with the same soil dataset, same TF, and subsequent spatial upscaling are shown in the second row (subplots 2a-c). Larger differences occur in regions with considerable subgrid heterogeneity in soil texture, where the hydraulic parameters associated with the mean textural information of the dominant class do not represent underlying variability. Absolute differences of more than 5% for both parameters occur in Florida, parts of Nebraska, and parts of the Southwest. The average 400 difference is -1.6 and -1.3 % for both parameter fields, and a greater variability can be found for θ s . The difference for the long-term ET flux is -3 % on average, with pronounced low values of less than -10 % in the aforementioned regions. The lower ET fluxes are due to the accumulated effects of lower porosity and lower hydraulic conductivity, which leads to decreased storage capacity and capability to meet evaporation water demands. It is important to keep in mind that these changes stem solely  (2006), and Vereecken et al. (1989Vereecken et al. ( , 1990, respectively. from using non-classified continuous soil textural information and aggregation of subgrid parameters, in contrast to using the 405 dominant soil type. In other words, variation originates from different methods of handling sub-grid variability. In the next step, we replaced the TF in MPR with another continuous TF based on the linear regression of the predictors of sand content, clay content, and organic matter (Saxton and Rawls, 2006) (maps 3a-3c). This TF is also available within reported range of values (0.01 to 6.6 % for organic matter content and 1.04 to 1.23 g/cm −3 for bulk density), upon which the regression for the TF was computed (Vereecken et al., 1989) is based on a small dataset of soil samples from Belgium. This does not cover the range of the values present in the chosen soil dataset (i.e., SG). The applicability of a given TF needs to be evaluated in the context of the utilized soil database (Wösten et al., 2001).

425
The model bias for a standard Noah-MP configuration in comparison to the reference product FLUXNET has been assessed in a previous study (Ma et al., 2017). They observed a spatial pattern in the long-term bias of ET in regards to prescribed leaf area index (LAI), which showed similarities to the pattern seen in Fig is that they directly investigated grid cell-specific parameters. As a result, the dimensionality of the parameter space linearly scales with the area of the model domain. Thus, only a few catchments using a spatially constant scale factor could be used.
Using MPR, the number of TF parameters remains independent of the size of the model domain, which allowed us to conduct 435 a more spatially comprehensive sensitivity analysis. At the same time, MPR requires TFs for every parameter, which indicates that the uncertainty due to the choice of TF cannot be neglected when conducting a sensitivity analysis.

Application with land-surface model HTESSEL
We used the land surface model HTESSEL (Balsamo et al., 2009) as an example to showcase the capability of MPR in coupling with state-of-the-art distributed environmental models. HTESSEL is the land-surface scheme used within the integrated 440 forecasting system developed at the European Center for Medium-Range Weather Forecasts (ECMWF). The model description and its configuration can be found in the Appendices B4, C4, the MPR configuration in the Supplement.
Similar to the application of Noah-MP (Fig. 5), we found differences in the long-term evapotranspiration (ET) flux of up to 15% over the Mississippi river basin (Fig. 6) when using different transfer functions to compute soil hydraulic properties.  model (Genuchten, 1980) of the hydraulic conductivity curve (Equation 1) and moisture retention curve (Equation 2): K s , α, n, l, θ r , and θ s . It is worth mentioning that this lookup table was originally derived for the FAO 2003 soil map (, CBL).
The most common values are 1.16e-6 ms −1 and 0.439 % for K s and θ s , respectively. These values correspond to the medium soil texture class in the lookup between soil hydraulic properties (i.e., hydraulic conductivity curves and moisture retention curve) and simulated ET at every location. Long-term ET values increases from 250 to 1450 mm per year along a gradient from the north-west to the south-east of the domain (Fig. 6 1c). We selected two sets of TFs from the literature to calculate soil hydraulic properties (Zacharias and Wessolek, 2007;Wösten et al., 2001). These TFs are applied to SG and can easily be implemented in MPR. After application of TF 2 (Zacharias and Wessolek, 2007), the parameter values θ s are reduced by about -9.0 % (Fig. 6 2b). The highest differences 460 occur again in Missouri and Kansas (around -30 to -47 % for θ s ). This TF does not contain Ks and the default Ks values are thus used. The change in θ s reduces long-term ET flux by about -5.0 % (Fig. 6 2c). This can be expected because soil moisture storage (i.e., porosity) is generally reduced (Fig. 6 2b). In turn, the amount of water available for plant transpiration is limited. The usage TF 3 (Wösten et al., 2001) for the estimation of K s and θ s reduces parameter values by approximately -6 % for Ks and -7 % for θ s on average over the entire domain. Applying TF 3 reduces porosity (Fig. 6 3b) in a similar order 465 of magnitude compared to TF 2. Additionally, saturated hydraulic conductivity is reduced at most by 16% resulting in reduced ET in comparison to the default setup (Fig. 6 3c). This reduction is not as strong as that of TF 2 because the reduced K s values increase the water holding capacity of the soil. highlights that the MPR-derived fields reduce the amplitude of the parameter values, but substantially increase the spatial variability. MPR-derived fields make more use of the spatial information of the soil dataset and lead to more realistic spatial parameter fields. It is worth mentioning that the spatial patterns for changes in ET do not correspond to neither changes in parameter fields nor the spatial pattern of the default ET values. This indicates the complex interplay between ET and soil hydraulic properties and calls for a deeper analysis of all MPR-derived soil parameters (i.e. also α, n, l and θ r etc.).

Differences between HTESSEL and Noah-MP over the Mississippi river basin
There are several differences between the simulations conducted with Noah-MP and HTESSEL that go beyond the fact that these are two are based on different mathematical models. First, the default simulations compute different amounts of longterm ET (compare 1c in figures 5 and 6). Both maps exhibit a similar spatial pattern, but the long-term ET flux for HTESSEL is approximately 20 % higher than that of Noah-MP. This might be due to the use of different forcing datasets NLDAS2 (Xia lookup-tables on the same dataset (SG) to rule out differences coming from different soil maps. While a decreased spatial variability especially for HTESSEL with only five active soil texture classes is found, the Cosby TF leads to a more consistent spatial pattern for Noah-MP. Additionally to the spatial pattern, the parameter values themselves are also different for both 490 models. θ s for Noah-MP varies around 0.46, and HTESSEL again shows slightly lower values of approximately 0.44. K s is on average around 2.5e-6 for Noah-MP and much lower with 1.16e-6 for HTESSEL. These striking differences are in agreement with Samaniego et al. (2017), where a more exhaustive model comparison was performed. This study again highlights the need for a common protocol to assess parameter uncertainty in distributed models.
Third, using other TFs than the default ones leads to reductions in long-term ET around less than 10% in magnitude for 495 both models (compare right column in figures 5 and 6). A similar magnitude of influence by varied soil parameters on ET has been reported previously by Livneh et al. (2015) for mHM in the Mississippi river basin. There is no consistent pattern between models in regard to where these changes manifest themselves. An example for that are the Nebraska Sandhills. While ET is generally increased by TFs in HTESSEL in this region, the opposite is the case in Noah-MP. Direct interpretation of the interplay of soil parameters in the soil water hydraulics is easier with Noah-MP due to the simpler mathematical model for soil 500 hydraulic parameterization. Noah-MP uses the Campbell parameterization to relate hydraulic conductivity to soil saturation (Campbell, 1974). In contrast, HTESSEL uses Mualem-van Genuchten parameterization (Genuchten, 1980), which leads to complex changes of moisture retention curves and hydraulic conductivity curves that are highly non-linearly impacted by In spite of demonstrated differences in model forcing, configuration and the process parameterization, MPR-derived parameter fields significantly changed long-term model output. The harmonization and reproducibility of parameter estimation across models through MPR opens up an avenue to a deeper understanding of the relationship between predictors, parameters and model parameterization.

Limitations and outlook
In spite of the versatility of MPR for its use cases, we acknowledge that limitations do exist, some of which are technical while others concern supported features. The limitations are: -The computational cost of MPR in large-scale applications and optimization applications is low due to a lack of parallelization. This will likely change in future releases by utilizing a hybrid MPI and Open-MP parallelization.

515
-MPR cannot yet handle transformation of different projected or geographic coordinate systems.
-Additionally, so far MPR only accepts the netCDF format as input or output. However, if the user community requests an extension to include GRIB-based, polygon-based, or HDF5-based input data, it may be implemented in future versions.
standards. MPR was tested with compiler GNU gfortran versions > 7.3 the NAG Fortran Compiler version 6.2, and Intel ifort version > 18.
We will invest further efforts in developing MPR so that scalability on high-performance computers (HPCs) and paralellization is improved. The need for support of advanced and massively parallel regridding and interpolation capabilities will likely lead to the integration of one of the existing remapping libraries used in GCM couplers in an upcoming version of MPR. The 525 coupling of MPR to the LSM HTESSEL by extraction of hard-coded and hidden parameters and the development of TFs in an ongoing project will likely serve as a template that can be adapted for other models.

Conclusions
Parameter regionalization enables the creation of seamless parameter fields for complex distributed models that can otherwise only be inferred through calibration or by default values, often obtained at inappropriate scales . MPR is 530 a framework that regionalizes parameters through the application of transfer functions and aggregations to any spatiotemporal coordinate system. In this study, we introduced a complete rewrite of the MPR framework (Samaniego et al., 2010;Kumar et al., 2013b) to overcome the limitations of previous implementations and comparable software (for example , Mizukami et al., 2017). MPR is able to introduce new flexibility to mHM and other models accepting distributed parameters through the support of multiple grid structures and a more flexible configuration of the parameter estimation process. It is capable 535 of reproducing effective parameter fields (Tóth et al., 2017) by applying TFs (Tóth et al., 2015), while also being able to remap/upscale the parameters onto every modeling unit (rectangular grids, HRUs, arbitrary shapes) required by the model. We demonstrated that parameter estimation not only exerts a strong influence on effective model parameter fields but also results in modified evapotranspiration simulated by land-surface models, even when MPR is applied to only two sensitive parameters.
The same holds true for models that use a tiling approach for handling subgrid heterogeneity.

540
The superiority of the MPR approach toward standard parameter estimation approaches was first demonstrated by Samaniego et al. (2010);Kumar et al. (2013b); Samaniego et al. (2017) and is now available for use with many other models. This is possible because MPR is designed to flexible, modular, and as easy to use. We provide an API that users can easily modify and that can be successfully tightly coupled to the hydrologic model mHM. As such, we invite implementation of further  (Tóth et al., 2017) at a given 250m and 1km resolution.

B3 Noah-MP
Noah-MP simulates the terrestrial water, energy, and carbon budget and estimates fluxes between various storage components 585 in the biosphere, lithosphere, and hydrosphere. Its predecessor Noah (Ek et al., 2003) was superseded by Noah-MP by implementing multiparameterization options and improved physics for various ecohydrological processes. For each grid cell, the vertical model structure was discretized into one canopy layer using a semi-tile approach, three snow layers, four soil layers, and an unconfined subsurface layer. mHM code was refactored and adapted to allow for the passing of global parameters from mHM to MPR, and to allow the effective parameter fields to be received.

C2 SoilHydroGrids
It is based on the SG dataset  at 250m and an aggregated 1km resolution. Linear functions as well as decision tree-based functions were used for the TFs (Tóth et al., 2015). A collection of relevant soil hydraulic parameters for 605 the European domain are provided in SHG. We selected the subdomain of the Netherlands for the parameter saturated hydraulic conductivity, as it shows a large degree of variability in this region. The MPR configuration file is attached in the Supplement.

C3 Noah-MP
We use the default WRF-Hydro parameterization process (version 3.6), at NCAR (2020), except for the radiative transfer option where a modified two-stream option was used. Meteorological model forcings were taken from the NLDAS-2 dataset (Xia 610 and NCEP/EMC, 2009). The 1/8 • spatial resolution and hourly temporal resolution of the forcing variables (air temperature, precipitation, specific humidity, wind speed, surface pressure, downward shortwave radiation, and downward longwave radiation) also constitute the chosen model resolution for Noah-MP. The SG dataset  was used to estimate soil parameters, whereas the MODIS-IGBP (FRIEDL et al., 2010) dataset was used to derive vegetation parameters. In the default setup, soil textural data was averaged over the soil column and classified into 12 classes using the Staff (1993) scheme. Finally,

615
the dominant type within a model grid cell was used. The exact parameters were then derived from class-specific default values provided in the lookup tables. These default values were derived by applying a set of TFs (Cosby et al., 1984) to the mean textural properties of the respective soil class. Vegetation parameters were also estimated based on the default values for each effective land cover class. The Noah-MP model version 3.9 was used and slightly modified to explicitly allow for two specific spatially distributed model parameters to be read. We used the default soil layering of horizon boundaries at depths of 0.1, 0.3, 620 0.6, and 1 m.
In addition to this default setup, we used MPR to estimate the soil parameters SATDK (soil saturated hydraulic conductivity) and MAXSMC (maximum soil moisture content). Noah-MP shows a substantial sensitivity to both of these parameters along a gradient of hydro-climate conditions in CONUS (Cuntz et al., 2016). The two parameters were estimated directly on a highresolution soil dataset using the following continuous TFs: TF 1 (Cosby et al., 1984), TF 2 (Saxton and Rawls, 2006), and 625 TF 3 (Vereecken et al., 1989(Vereecken et al., , 1990. TF 1 is used in Noah-MP as a default option, TF 2 was later introduced as an option in Noah-MP version 4.0, and TF 3 was chosen in this study to demonstrate the effect of a TF based on soil samples from outside the study domain. The arithmetic mean was used to upscale the parameters to the model resolution, except for the vertical scaling of K s along the soil horizons for which the harmonic mean was used. The parameters REFSMC (soil moisture content at field capacity) and WLTSMC (soil moisture content at wilting point) were rescaled by the ratio of the default and modified

635
Meteorological forcing data were taken from the ERA5 dataset (ECMWF, 2019). The 1/4 • spatial resolution and 3-hourly temporal resolution of the forcing variables (air temperature, precipitation (rain and snow), specific humidity, wind speed, surface pressure, downward shortwave radiation, and downward longwave radiation) also constitute the chosen model resolution for HTESSEL. 1 abstract interface function t ra n s f e r _ f u n c _ a l i a s (x , param ) result ( func_result ) ! import the double precision kind specification and custom type import dp , I n pu t F i e l d C on t a i n e r ! > an array containing the predictor variables ( access values through ' data_p ' property ) 6 type ( I n p u t F i e l d C o n ta i n e r ) , intent ( in ) :: x (:) ! > an array containing the TF parameters real ( dp ) , intent ( in ) :: param (:) ! > the resulting TF result real ( dp ) , allocatable :: func_result (:) 11 ! ! allocate the func_result to the size of the predictors ( all have the same size ) ! allocate ( func_result ( size ( x (1) % data_p ) ) ) ! ! enter the TF function here ! func_result = x (1) % data_p + x (2) % data_p * param (1) 16 end function t r a n s f e r _ fu n c _ a l i a s end interface Figure D1. Template for user-defined TFs as an abstract interface.
We used a process parameterization based on the default configuration presented in the development branch CY47R1 of 640 HTESSEL (nemk_CY47R1.0_v6b_cmflood_mpr). In addition to this default setup, we used MPR to estimate the six soil parameters of the Mualem-van Genuchten model for the hydraulic conductivity curve (Equation 1) and soil moisture retention curve (Equation 2). These parameters were estimated directly on the SG dataset  using the following TFs: categorical TF 1 based on lookup-table values (HYPRES, 1997), continuous TF 2 (Zacharias and Wessolek, 2007) with the estimation of only the four parameters of Equation 2, and continuous TF 3 (Wösten et al., 2001). Soil textural data was 645 averaged over the soil column for TF 1, classified into 7 classes using the HYPRES soil texture triangle (HYPRES, 1997;Wösten et al., 1999), with some additions for organic soils. Finally, the dominant type within a model grid cell was used. The exact parameters were then derived from class-specific default values provided in the lookup tables. TF 1 is used in HTESSEL as a default option. The arithmetic mean was used to upscale the parameters to the model resolution, except for the vertical scaling of K s along the soil horizons, for which the harmonic mean was used.

650
The HTESSEL model version CY47R1 was used and modified to explicitly allow for the reading of spatially distributed model parameters. We used the default soil layering of horizon boundaries at depths of 0.07 m, 0.28 m, 0.1 m, and 2 m. Due to limitations in the HTESSEL solver for the soil physics processes, we enforced vertically homogeneous soil parameters.
The HTESSEL model was run for 8 years at a daily time step from 1979 to 1986. We allowed the model to run for an entire period and used the resulting state variables as initial conditions for the final run. The final evaluation period covered the years 655 1 abstract interface real ( dp ) function u ps ca l e _ f u n c_ a li as ( array , weights , p ) ! import the double precision kind specification import dp ! > the array of subgrid values ( no missing values ) 6 real ( dp ) , dimension (:) , intent ( in ) :: array ! > the array of weights ( same shape as array ) real ( dp ) , dimension (:) , intent ( in ) , optional :: weights ! > an optional parameter passed to the function ( power mean ) real ( dp ) , intent ( in ) , optional :: p   x of the array with size n (with indices i) are passed to the operator with weights w of the same size.