Physical Constraints for Improved Soil Hydraulic Parameter Estimation by Pedotransfer Functions

Global land surface models use spatially distributed soil information for the parameterization of soil hydraulic properties (SHP). Parameters of measured SHP are correlated with easy‐to‐measure soil properties to construct general pedotransfer functions (PTFs) used to predict SHP from spatial soil information. Global PTFs are based on a limited number of samples yielding highly variable and poorly constrained SHP. The study implements a physical constraint, soil‐specific capillary length, to reduce unphysical combinations of SHP. The procedure fits concurrently soil water retention and capillary length using the same parameters. Results suggest that meeting the capillary length constraint has minor effects on the goodness of fit to soil water retention data. Constrained SHP were applied to represent 4 years of lysimeter fluxes yielding evapotranspiration values in close agreement with measurements relative to slight overestimation by unconstrained SHP. The procedure was applied for testing constraint SHP at a regional scale in New Zealand using the surface evaporation capacitance model and Noah‐MP for detailed simulations of land surface processes. The use of constrained SHP in both models yields higher surface runoff in agreement with observations (unconstrained SHP severely underestimated runoff generation). The concept of constrained SHP could be extended to include other physical constraints to improve PTFs, for example, by consideration of vegetation cover and soil structure effects on infiltration.


Introduction
The quantification of near-surface hydrologic processes requires information on spatially distributed parameters characterizing the soil hydraulic properties (SHP). These parameters are inferred from fitting models to measured values of soil water characteristics (SWC) that relate volumetric soil water content θ with the capillary head h. The same SWC parameters are also used to describe the unsaturated hydraulic conductivity function K(θ) often anchored by a measured value of saturated hydraulic conductivity for the same soil type or textural class. The parameter values for the SHP are estimated either (i) by minimizing deviations between measured and parametrized water content capillary head (or conductivity) values or (ii) by using pedotransfer functions (PTFs) that link easy-to-measure soil properties (textural fractions, bulk density, organic carbon, or pH) with the hydraulic properties. The resulting SHP are not only important for describing the SWC for a location but are also used to determine other physical properties well beyond the information considered in the estimation process. For example, estimating SWC requires determination of parameter values that describe the pore size distribution and a characteristic pore size; these values play an important role in other physical processes (thermal conductivity, water infiltration, and capillary flow towards an evaporating surface).
PTF-estimated SHP often exhibit large scatter; Gutmann and Small (2007) have shown that soil textural classes explain only 5% of the variance of SHP (in other words, SHP variations exceed expected differences due to soil texture classification). To overcome ambiguities in the determination of SHP based on goodness of fit to SWC only, we report a method for injecting additional physical constraints to reduce occurrence of unphysical parameter combinations and thus potentially improve land surface representation. The primary physical constraint used here is a soil-specific capillary length that governs capillary flow to an evaporating surface (Lehmann et al., 2008). Experimental evidence based on evaporation dynamics from drying soils (Gardner & Hillel, 1962;Jackson et al., 1976;Lehmann et al., 2008;Salvucci, 1997) and evaporation rates from shallow water tables (Gardner & Fireman, 1958;Hellwig, 1973;Mengistu et al., 2018;Rasheed et al., 1989;Shokri & Salvucci, 2011) shows that the maximum depth that sustains capillary flow from wet soil to an evaporation surface rarely exceeds 1 m. This is depicted in Figure 1a, based on experimentally determined maximum length of capillary flow paths. The evaporative capillary length is characterized by the same SWC parameters for a given soil (see Lehmann et al. (2008) and analyses in section 2).
As shown in Figure 1b, for a subset of the experiments where SHP are reported, the predictions based on SWC parameters estimates are in good agreement with the inferred length of capillary flow paths as explained shortly. The hypothesis is thus that injecting this physical constraint would improve the estimation of the parameters shared by the SWC and the characteristic length for capillary flow at a given soil.
The specific objectives of this study were to (1) illustrate how the standard and unconstrained SHP (USHP) estimation may result in unphysical parameter combination (relative to the expected characteristic capillary length for a soil type), to (2) propose a parameter estimation procedure that links SWC and characteristic length information, and (3) to evaluate the potential improvement of constrained SHP (CSHP) for near-surface process representation. The concept of evaporation characteristic length (Lehmann et al., 2008) serves to highlight some of the shortcomings of present PTF-based SHP and offers a means for rectifying unphysical parameter combinations. In section 2 we introduce the parametric representation of the characteristic length for capillary flow. In section 3, we present the data and procedure for constraining the SHP. The results section highlights applications and tests of the procedure using lysimeter measurements; we then discuss potential applications for larger scales relevant to land surface models (LSMs) using the water balance of New Zealand southern island, followed by a summary of the results and conclusions.

Evaporation Characteristic Length: A New Physical Constraint
Evaporation dynamics from drying soils are often characterized by an initially high and relatively constant evaporation rate E 0 (Stage I evaporation) supported by capillary flow to the surface (Prat, 2002;van Brakel, 1980), followed by a decreasing rate (Stage II evaporation). In Stage II, the vaporization plane (where phase change occurs) recedes into the porous medium and the evaporation rate is controlled by diffusion from the vaporization plane through the porous medium. During Stage I evaporation, the atmospheric demand for phase change at the surface is supplied by capillary flow from within the wet porous medium (Scherer, 1990). Lehmann et al. (2008) formulated an intrinsic quantity that marks the depth over which capillary flow is maintained and supplies the evaporating surface. This length L is deduced from the SWC parameters (the difference between air entry value h b and the pressure at critical capillary head h crit ) and reflects the conditions for interruption of capillary flow paths. Lehmann et al. (2008) determined the . The shaded area marks the range associated with experimental results. The two continuous black lines describe theoretical predictions using different soil hydraulic parameter data sets (see section 3.2). (b) For a subset of experiments with reported soil hydraulic properties (see references in legend and in text), the characteristic length for capillary flow L C was estimated using the approach derived in Lehmann et al. (2008). maximum extent of this capillary flow path for a given soil by defining a characteristic length L C of evaporation (Lehmann et al., 2018): with Stage I evaporation rate E 0 and the unsaturated hydraulic conductivity for a critical capillary pressure h crit that marks the disruption of capillary flow paths. A key point in the derivation in Lehmann et al. (2008) is the use of the SWC information (or the van Genuchten parameters) to define this length. These are the parameters n, α, and m = 1 − 1/n that appear in equation 1. For completeness, the hydraulic functions linking water content θ, capillary pressure h, and hydraulic conductivity K based on the model of Mualem (1976) and van Genuchten (1980) are given in equations 2a and 2b: with effective water content Θ, residual θ res and saturated water content θ sat , the saturated hydraulic conductivity K sat , and the tortuosity factor τ set to 0.5. As shown in Figure 1a and in supporting information S1, (i) L C is large for intermediate textures and small for coarse and fine textures (due to limited capillary driving forces for coarse media and high viscous resistance in fine textures, respectively) and (ii) L C rarely exceeds 1 m for evaporation rates E 0 of several millimeters per day (see experimental results in Figure 1).
It was shown in previous studies (see Figure 1b) that for well-defined parameter values of n, α, and K sat , equation 1 provides a reliable estimate of the maximum length L C of capillary flow. However, for certain parameter combinations, unrealistic high values of several meters are computed using equation 1. To differentiate between "permissible" and "forbidden" parameter combinations and lengths L C , we use average properties of thousands of soil measurements. These measurements provide the range for realistic parameter combinations. To quantify these parameter combinations, we define α and K sat as function of n (for definition of relationship α(n) and K sat (n) see equations in supporting information S2). We chose n as independent variable because it increases systematically with soil "coarseness" that makes n a useful representative or surrogate for soil texture (n close to 1 for clay and n > 2 for sand). To obtain estimates of characteristic length L C that are within the limits defined by experimental results, we insert the empirical relationships α(n) and K sat (n) into equation 1 and obtain the analytical expression for L C denoted as target value of characteristic length L t for a given soil type or textural class (characterized by n): with parameters a, b, c, and d used to describe the correlation between n, α, and K sat and with a minimum value of parameter n min (see details in supporting information S2). The requirement for fitting SWC is that the characteristic length L C calculated in equation 1 should also match the analytical expression for the same soil texture given in equation 3 to fulfill the physical constraint of the characteristic length of evaporation. How this is done is explained next.

Constrained Estimates of SWC Parameter Values
For the joint parameter estimation, we define an objective function f(α, n) that provides a weighted sum of the squared residuals between estimated and measured SWC water content and between L C estimated based on equation 1 and the soil type target value given in equation 3: with the number N of pairs of the measured SWC, the estimated water content θ for a certain capillary head value h i , and the corresponding target (i.e., measured) value θ t , the calculated characteristic length L C and the target value L t , and the weights ω θ and ω LC assigned to the residuals of the water content and the characteristic length for capillary flow, respectively. To solve the nonlinear least squares problem (minimizing equation 4), we applied the Gauss-Newton method (see supporting information S3).

SWC Experimental Data
Applying physically constrained fitting of SHP as expressed in equation 4, we expected that the fitted values of n and α would differ from values for the unconstrained fitting. We have used the UNSODA database (Nemes et al., 2001) selecting reported values with at least six SWC data points (measured in the lab, not in the field) and with information of saturated hydraulic conductivity. For the resulting 377 SWC data sets, we estimated the parameters n, α, and L C using minimum and maximum water content as estimate for residual and saturated water content (θ res and θ sat ). In average 11 pairs of water content θ and capillary head h were available to estimate the parameter values. For 368 of the data sets, the soil textural class was also reported and enabled association with mean values for 12 soil textural classes (note that for silt, sandy clay, and clay loam, less than 10 data series could be used). The average soil textural properties could also be compared with the mean values reported in Carsel and Parrish (1988; they used samples collected for soil survey data reports) and for all UNSODA data estimated by Leij et al. (1999) using neural network predictions (Schaap & Leij, 1998) and other PTFs (Rawls & Brakensiek, 1985).

Global Maps of SHP
To simulate processes in the vadose zone at the global scale, locally distributed SHP are required. The SHP are often estimated using PTFs based on basic textural information as provided in the SoilGrids database (Hengl et al., 2017). We apply here the Rosetta3 PTF (Zhang & Schaap, 2017), estimating n, α, and K sat , from bulk density, sand, silt, and clay content (for topsoil with depth 0.0 to 0.3 m) using a neural network. The global maps with resolution of 0.1°are shown in Figure 2. The characteristic length was calculated with equation 1 using a potential evaporation rate E 0 of 4 mm/day. This value was close to the global average over the terrestrial surfaces (3.8 mm/day) as computed with the approach of Jensen and Haise (1963) using incoming shortwave radiation data from CERES EBAF (https://ceres.larc.nasa.gov/) and temperature data provided by the Earth System Research Laboratory (https://www.esrl.noaa.gov/).

Simulation of Lysimeter Flux Data
An important question in the core of the study is whether CSHP would improve model performance. The systematic evaluation of this question would require dedicated studies and tests within spatially distributed systems (see section 3.5). To provide a proof of concept, we opted for testing the constrained SWC fitting on lysimeter measured soil water dynamics by comparing simulations using CSHP and USHP with observations. The lysimeter data were obtained from a site from the Tereno network (Pütz et al., 2016) called "Wüstebach" in Germany for the period January 2012 until end of April 2016 (~1,500 days). The lysimeter was installed in a clearing in the forest at elevation of 610 m. Samples from vicinity of the lysimeter were taken to measure SWC and saturated hydraulic conductivity. Evapotranspiration ET was delineated from mass balance, using measured lysimeter mass and rainfall. Potential evapotranspiration ET 0 values were estimated using Penman-Monteith based on data recorded by a nearby weather station. All data were provided by H. Bogena (see acknowledgments section).
Evaporation and other processes of the water cycle were simulated using the surface evaporation capacitor (SEC) introduced in Or and Lehmann (2019). The SEC represents a soil layer with thickness equal to characteristic length L C , considering L C as maximum depth from where water can be extracted by evaporation (water percolating to deeper soil depth is protected from evaporation). The SEC concept considers daily dynamics of Stage I and Stage II evaporation E and redistribution Q to deeper soil layers, when water content in the capacitor is larger than the critical water content θ crit (expression for θ crit is given in supporting information S1). For vegetation-covered surfaces as in the case of the Wüstebach-lysimeter, evaporation of intercepted water is considered as well. As a simple estimate of transpiration T, we assume that water percolating from the bottom of the capacitor layer Q is converted to transpiration T = min(Q, ET 0 − E) and deep percolation D = Q − T (i.e., Q is converted entirely to transpiration T when Q is smaller than the difference between potential evapotranspiration ET 0 and surface evaporation E).
As shown in equation 1, the characteristic length is a function of n and α. Differences in L C due to constrained fitting will affect the storage dynamics as follows: the shorter L C , the higher is the change in water content after rainfall P, assuming that the rainfall is distributed homogeneously within the capacitor. Higher water content corresponds to higher discharge and more Stage I evaporation.

Large-Scale Water Balance Representation (New Zealand)
SHP applied in LSMs may result in overestimation of capillary pull (i.e., overestimation of L C ). To quantify the effect of constrained and unconstrained parameter values, we calculated the water balance for the region of New Zealand for several years. We selected the southern island of New Zealand due to availability of hydrologic information for a small area and many sites with a characteristic length of several meters (see section 4). We took as input the SHP obtained with Rosetta3 and recalculated the estimated SHP constrained by L C .
In a first step we calculated the water balance for the years 2001-2016 using the one-dimensional surface evaporation capacitance (SEC) concept (Or & Lehmann, 2019). The SEC requires input information on potential evapotranspiration, rainfall, and vegetation. The leaf area index is used for estimation of canopy rainfall interception and shading effects of canopy on potential evaporation rate. For leaf area index we have used data from the Copernicus Global Land Service (https://land.copernicus.eu/global/products/lai). Rainfall data were taken from MSWEP database (Beck et al., 2019), and monthly values of mean potential evaporation (climatic average) were estimated by the method of Jensen and Haise (1963).
In a second step, to check if constrained parameterization of SHP (in the following called CSHP in contrast to unconstrained parameterization USHP) affects the dynamics of a complex LSM, we used the LSM Noah-MP (Niu et al., 2011) and applied it for the years 2004 and 2005. Noah-MP describes mass and energy fluxes in interacting soil and atmosphere including root zone water uptake and canopy resistance. The soil water dynamics including infiltration excess were calculated for four soil layers by solving the Richards equation (Richards, 1931). By default, Noah-MP uses hydraulic functions according to Brooks and Corey with parameter values assigned to soil textural classes by Cosby et al. (1984). Here the model was adapted to use (i) van

Evidence of Unphysical Values in Present (Unconstrained) Global SHP Databases
The hydraulic parameters of terrestrial surfaces with a resolution of 0.1°(≈10 km) were estimated with Rosetta3 PTF. The results are shown as maps in Figure 2.
Inserting values for n, α, and K sat in equation 1, a map of the characteristic length L C for capillary flow was generated. In Figure 2, two features should be highlighted: (i) The sandy desert region of Africa is characterized by high values of the hydraulic parameters, and (ii) the soils of northern regions are mainly loam and sandy loam with a high characteristic length. As was mentioned in section 2 and supporting information S2, the characteristic length is short in coarse (narrow pore size distribution) and fine textures (high flow resistance) but large in intermediate texture as manifested in Figure 2d.
In northern latitudes, the characteristic capillary length deduced from equation 1 exceeds several meters and is thus considered unphysical (see Figure 1 above). For about 30% of the global terrestrial surfaces, we estimate L C values that are larger than 1 m with a maximum value of 8.8 m. In Figure 3a, the distribution of the characteristic length L C is shown as function of hydraulic parameter n and reveals that values were larger than expected based on previous studies (see Figure 1a). Interestingly, SHP values for regions of the African desert (Figure 3b) follow the expected trend with short values for coarse and fine and higher L C values for intermediate textures.
As illustrated in Figure 3a, L C values estimated from mean SHP of all UNSODA samples (Leij et al., 1999) and for Carsel and Parrish (1988) parameters are lower than 1 m (with the exception of the silt value for UNSODA data that was estimated based on three SWC samples only). To relate the observed large values at global scale to "expected" values calculated from mean values of the 12 soil textural classes, we describe the relationship between average values of n, α, and K sat , using equations given in supporting information S2 to calculate the characteristic length as function of shape parameter n given in equation 3.
A potential explanation of the low α values in UNSODA database (and the corresponding high values of the characteristic length L C ) is the potential role of soil structure not explicitly accounted for. SHP are sensitive to the presence of soil structure (biopores and aggregates) for the same soil texture (Or, 2019). A partial support for this conjecture is provided by the SHP for African desert soils with limited vegetation and limited structure forming biological activity in which the calculated capillary characteristic length is in agreement with theoretical values. Conversely, poorly CSHP and overestimation of the characteristic length for soils with structure-forming processes suggest a role for this ubiquitous factor (Weynants et al., 2009) not presently considered in PTFs (Gutmann & Small, 2007;Or, 2019).

Fitting SWC With Capillary Length as a Constraint
The set of 377 measured SWC curves in UNSODA database was fitted with weights NÁω θ þ ω LC ¼ 1 assigned to the residuals of the water content (with N the number of pairs of water content and capillary head for a certain sample) and the characteristic length, respectively. We tested different weigh ratios ϖ ¼ ω LC =ω θ ranging from 0 (no physical constraints) to 0.02. The small weight ratio can be explained with the following example: For 10 sample points with deviation between measured and fitted water content θ of 0.01 m 3 /m 3 (1%) and a weight ω θ of about 1/10, the total residual (sum of 10 quadratic differences) is 0.0001; when the difference between calculated L C and the target expression L t from equation 3 is about 1 m, the quadratic deviation in L C is 10,000-fold the deviation or one θ value. By assigning a weight ω LC ¼ 0:001ω θ to the deviation in L C , both errors would be weighted equally strong.
As shown in the supporting information S3, for ϖ values larger than 0.004 the error in fitted water contents was doubled compared to the unconstrained fitting and was considered as too high. In the following, we chose a weight ϖ equal to 0.002. One illustrative example on the effect of physical constraint is shown in Figure 4a, with a characteristic length dropping from 1.2 to 0.4 m by considering the deviation of the characteristic length L C from target value L t . As shown in Figure 4b, without physical constraint many characteristic length values are larger than 1.5 m (43 of 377 samples). No characteristic length L C >1.5 m was obtained for the constrained fitting. The average L C for constrained fitting was 0.48 m compared to 0.84 m for the unconstrained fitting.
The effect of constrained fitting on values of n and α is not the same for all samples as shown in Figures 4c and 4d. In average, the geometric mean for n increases from 1.97 to 2.09 and α increases from 2.06 to 2.12 m −1 when SWC was fitted considering physical constraints. With respect to the accuracy of the fitting, under the physical constraints reduced the mean deviation of the capillary characteristic length from 0.63 to 0.19 m at the "cost" of increasing mean deviation of SWC in terms of water content from 0.030 to 0.045 m 3 /m 3 .

Effect of CSHP on Representation of Lysimeter Hydrologic Fluxes
The time series of lysimeter water balance information from Wüstebach for a period of more than 4 years was simulated using the SEC concept (Or & Lehmann, 2019). The parameters for the SWC were estimated with and without the characteristic length physical constraint (Figure 5a). The characteristic length deduced for the soil was 0.4 m for the constrained parameterization (CSHP) and 1.1 m for the unconstrained scenario (USHP). This difference resulted in shallower capacitor depth for the constrained case and thus affecting storage dynamics as shown in Figure 5b. Both SEC simulations (CSHP and USHP) were in good agreement with the measurements, with slightly better performance for the constrained parameters. The model overestimates the measured cumulative evaporation fluxes by 280 mm for USHP and 60 mm for CSHP. The difference in model-inferred evapotranspiration was thus 220 mm with higher evapotranspiration and lower discharge for USHP. The differences in simulated evapotranspiration fluxes stem from higher transpiration values for unconstrained parameters (see section 3.4 on the partitioning of rainfall water to surface evaporation, transpiration, and deep drainage). Due to larger water contents after rainfall in case of constrained fitting with shorter L C , more water Q will percolate to deeper soil layers and will be converted to deep drainage D and not to transpiration T. We thus conclude that the constrained parameterization performed slightly better than the original unconstrained parameters.

Discussion on Effect of CSHP in LSMs
As shown in section 4, the values of the characteristic length L C calculated from the Rosetta3 database often exceed the physical range of values for the loamy soil types (see Figures 1 and 3) for various regions around the globe. We expect that LSMs representation of near-surface processes based on these parameters would lead to overestimation of surface evaporation and potentially stronger coupling with shallow groundwater due to the unrealistically large capillary flow length. However, how could such a method be applied in a spatially distributed manner to update the presently used values of SHP? A potential correction could be implemented by considering the unconstrained parameter values for locations from Rosetta3 (or another global SHP data set) to deduce 10 values of the SWC (water content θ i = θ res +[0.05+(i − 1)/10](θ sat − θ res ) and the corresponding capillary pressure). These pairs of (θ i , h i ) (assumed i = 10 in the following example) are considered as local "measurements" for readjusting SWC with the constraint that the L C values are close to the expected values L t (equation 3) for the local soil type. We implemented such a local adjustment procedure for SHP from New Zealand as seen in Figure 6 where the original characteristic length values exhibited values in excess of 6 m. The refitting of the θ i − h i pairs subject to the proposed physical constraint L C (implementing a weight ratio ϖ ¼ ω LC =ω θ ¼ 0:5 for the adjustment with such high L C values) reduces the characteristic length values ( Figure 6a) and changes the spatial distribution of parameter values (Figures 6b and 6c). Evaluating the benefits of constraining SHP especially at large scales on near-surface process modeling is not a simple task. In the following, we test the resulting SHP and compare observations and simulations for the southern island of New Zealand as a test case. We distinguish between simulations with the capacitor model and a full LSM approach (Noah-MP).

Simulations of Water Balance With Surface Evaporation Capacitance (SEC) Concept
We calculated components of the surface water balance using the surface evaporation capacitance (SEC) concept (Or & Lehmann, 2019). The surface evaporation was very similar (228 mm for CSHP and 230 mm for USHP). More relevant are differences with respect to surface runoff generation due to saturation excess (R off ). Figure 7 depicts differences between unconstrained (U) and constrained (C) parameters for redistribution below the evaporation active capacitor (Q) and surface runoff generation. The large variability in rainfall distribution between the east and west sides of the Southern island of New Zealand is reflected by a highly variable spatial distribution of water balance that is not reflected by the mean values. For the simulations using CSHP, we observe runoff in the high rainfall regions, whereas no runoff is predicted for USHP simulations. The difference in runoff averaged over the entire New Zealand was relatively small (20 mm); however, locally, values in excess of 400 mm/year were estimated with constrained parameters in the heavy rainfall regions in the West Coast of the Southern Island.
To assess if the surface runoff on West Coast for simulations with CSHP were in agreement with observations, we considered runoff from the Haast River catchment (data source see section 3.5), dividing the measured discharge by the catchment area of 1026 km 2 (Jowitt, 1999). River discharge data were subdivided into direct runoff and base flow using the filtering approach of Eckhardt (2005). Observations suggest that direct runoff was about 1.8 m out of the cumulative 17 m discharge (for 3 years); this estimate of runoff is in qualitative agreement with 1.0 m for SEC-simulated runoff excess using the constrained parameters. In other words, we demonstrate that direct surface runoff occurs and can be predicted using CSHP but not with USHP (with "deeper" capacitor and lower values of the parameter α).

Simulations Using LSM (Noah-MP)
The findings of increased surface runoff for constrained parameters simulated with capacitor model for New Zealand were evaluated using the Noah-MP LSM. As shown in Figure 8, the simulations with constrained The river discharge measured at local station was then subdivided into direct runoff (shown in blue) and baseflow using the approach of Eckhardt (2005). The gaps stem from missing values in the precipitation data series. The saturation excess simulated with the surface evaporation capacitor (SEC) using constrained parameters is shown in red.

Water Resources Research
parameter values resulted in higher surface runoff compared to the unconstrained parameterization and were thus in agreement with the surface runoff deduced from measurements and with simulations using the SEC model. Quantitatively, the simulations with Noah-MP for a period of 2 years produced 455 mm for the constrained and 110 mm for unconstrained parameterization. Both values are smaller than 1100 mm deduced from the measured runoff time series. Note that the surface runoff was not directly measured but deduced from filtering the river discharge time series, so the "true" surface runoff value is unknown and it cannot be concluded on the accuracy of the model predictions. The soil evaporation was simulated as well and manifested large effects of the constrained parameterization with more evaporation using USHP (see supporting information S4).
The inset in Figure 8 illustrates that modeled surface runoff patterns by the Noah-MP and SEC models are generally similar; however, there are also differences between the two models. These differences are attributed to three different aspects: (i) differences in the spatial resolution of the simulations (1 km for Noah-MP, 20 km for SEC), (ii) different input data (using different rainfall data will have a considerable effect), and (iii) differences in the process representation (physics) implemented in the two models (Noah-MP considers several layers with soil evaporation occurring only from top 0.15 m, whereas the SEC considers the thickness of the soil layer contributing to soil evaporation as a function of soil texture). Despite these differences, both model approaches exhibit similar effects of constrained parameterization, such as higher surface runoff for physically constrained characteristic length.
The sensitivity of Noah-MP to input hydraulic parameters for simulations of surface runoff has been shown in previous studies (Cuntz et al., 2016;Mendoza et al., 2015). In a detailed and systematic sensitivity analysis by Cuntz et al. (2016), it has been shown that not only the soil saturated hydraulic conductivity and porosity impart sensitivity for various hydrologic processes but also the shape parameter of the SWC curve (see parameter "BB" in Figure 4 of Cuntz et al., 2016). We thus expect that the effect of physical constraints on the resulting shape parameters of SWC (see Figure 4 above) plays a role in the simulation results using Noah-MP (and potentially other LSMs). For the analysis of southern island, the surface runoff at the West Coast with unconstrained (a) is much smaller than for the constrained parameter values (b). Note that the coarse resolution of climatic data (1/4°) is manifested in square-shaped structures in the runoff maps. (c) Also for the Haast catchment, the surface runoff for constrained parameters (red line) is higher than for unconstrained parameters (blue) and in better agreement with the values deduced from measurements by filtering (black line).

Summary and Conclusions
LSMs rely on estimates of hydraulic properties that are needed to represent near-surface water and energy fluxes. However, these parameters that are often deduced from soil textural properties and PTFs are poorly constrained and highly variable. Among the strategies for improving SHP for Earth systems and LSM applications, we propose the use of physical constraints in addition to a few SWC data sets. An example of such physical constraint is the evaporation characteristic length L C for capillary flow that rarely exceeds 1 m while combinations of USHP predict much larger L C values in some cases. The study illustrates that it is feasible to simultaneously honor measured SWC and L C with minor loss of accuracy of fitting SWC. The results can be summarized as follows: • The deduced capillary characteristic length based on present SHP (Rosetta3) is overestimated (with L C value greater than 1 m) for about 30% of global terrestrial area. • We capitalized on correlations between shape parameters n, α, and saturated hydraulic conductivity K sat to derive an analytical function of L C (used as target value L t ) for implementing the physical constraint for each soil type. • For 377 soil samples with reported SWC data, the procedure of refitting SHP under the physical constraints reduced the mean deviation of the capillary characteristic length from 0.63 to 0.19 m at the "cost" of increasing mean deviation of SWC in terms of water content from 0.030 to 0.045 m 3 /m 3 . • Fitting SWC with the proposed physical constraint resulted in a reduction of the mean value of capillary characteristic length from 0.84 to 0.48 m. • Using physically constrained fitting of the SHP for a lysimeter study reduced the SEC-estimated evapotranspiration (cumulative over more than 4 years) by 220 mm relative to USHP and is in better agreement with measurements. • Preliminary implementation of the physical constraints at a regional scale (New Zealand) shows significant differences in SHP parameters between constrained and unconstrained fitting; models using the CSHP predicted much more surface runoff than using unconstrained parameterization in agreement with catchment observations.
The concept and proposed procedure of constraining SHP estimation from PTFs could be extended to include additional physical constraints, for example, the injection of dynamic parameters such as the infiltration-based "time of ponding" (the time at which soil infiltration capacity for a rainfall rate is exceeded and surface ponding ensues). Such a physical constraint is of particular interest due to its sensitivity to the presence of soil structure known to be absent in most PTFs (Gutmann & Small, 2007;Or, 2019). The demonstrated addition of capillary characteristic length to the estimation of SHP and additional potential physical constraints that link soil textural class and soil cover (vegetation and land use) may offer a means for improving the quality of pedotransfer-derived SHP used in land surface, Earth systems and large-scale hydrological models.