On the Uncertainty Induced by Pedotransfer Functions in Terrestrial Biosphere Modeling

Hydrological, ecohydrological, and terrestrial biosphere models depend on pedotransfer functions for computing soil hydraulic parameters based on easily measurable variables, such as soil textural and physical properties. Several pedotransfer functions have been derived in the last few decades, providing divergent estimates of soil hydraulic parameters. In this study, we quantify how uncertainties embedded in using different pedotransfer functions propagate to ecosystem dynamics, including simulated hydrological fluxes and vegetation response to water availability. Using a state‐of‐the‐art ecohydrological model applied at 79 sites worldwide, we show that uncertainties related to pedotransfer functions can affect both hydrological and vegetation dynamics. Uncertainties in evapotranspiration, plant productivity, and vegetation structure, quantified as leaf area, are in the order of ∼10% at annual time scales. Runoff and groundwater recharge uncertainties are one order of magnitude larger. All uncertainties are largely amplified when small‐scale topography is taken into account in a distributed domain, especially for water‐limited ecosystems with low permeability soils. Overall, pedotransfer function related uncertainties for a given soil type are higher than uncertainties across soil types in both hydrological and ecosystem dynamics. The magnitude of uncertainties is climate‐dependent but not soil type‐dependent. Evapotranspiration, vegetation structure, and plant productivity uncertainties are higher in water‐limited semiarid climates, whereas groundwater recharge uncertainties are higher in climates where potential evapotranspiration is comparable to precipitation.

. A major source of uncertainty when solving the Richards equation for the variably saturated critical zone (the region where vegetation typically accesses water) is the highly nonlinear relationship between soil water content and water potential (i.e., the soil water retention curve) and the dependence of hydraulic conductivity on the degree of saturation (i.e., the hydraulic conductivity function) (Assouline & Or, 2013). Models typically adopt the Brooks and Corey (BC) or the Mualem-Van Genuchten (MvG) formulations (Brooks & Corey, 1964;Mualem, 1976;Van Genuchten, 1980) to describe these nonlinear relationships. The estimation of model parameters for the BC and MvG models requires extensive data collection and is commonly performed, in terrestrial biosphere models, by means of pedotransfer functions (PTFs) (e.g., Assouline & Or, 2013;Patil & Singh, 2016;Van Looy et al., 2017;Vereecken et al., 2010). PTFs are essentially empirical models that use easy-to-measure soil properties (mostly soil texture, bulk density, and organic content) to predict the saturated hydraulic conductivity and the parameters of either the BC or the MvG formulations. PTFs use different statistical techniques, from multiple linear regression to advanced machine learning, to obtain those parameters based on a generally limited set of observations (e.g., Patil & Singh, 2016;Van Looy et al., 2017;Wösten et al., 2001;Y. Zhang & Schaap, 2019). Some soil textures (e.g., clays) are typically under-sampled and observations are usually obtained within a limited geographic region, leading to highly variable predictions across PTFs trained with different data sets (e.g., Wagner et al., 2001) and limited transferability across geographic locations (Gupta et al., 2020;Vereecken et al., 2016). The use of PTFs in both land surface and terrestrial biosphere modeling is ubiquitous (e.g., Cooper et al., 2020;Van Looy et al., 2017;Vereecken et al., 2016), and thus understanding the implications of their use is crucial. Global scale application of such models, in particular, rely heavily on PTFs (e.g., Dai et al., 2019), as specific soil hydraulic properties are available only locally.
The disagreement across PTF predictions is expected to impact hydrological dynamics, and this has been confirmed by a number of studies, where soil moisture (e.g., Deng et al., 2009;Loosvelt et al., 2011;Pinnington et al., 2021), evapotranspiration (e.g., Weihermüller et al., 2021), runoff production, and flood generation have been found to be strongly affected by soil hydraulic parameters (e.g., Mohajerani et al., 2021;W. Sun et al., 2016). The uncertainty introduced by PTFs in the hydrological dynamics has also been found to be modulated by small scale topography (e.g., Mohajerani et al., 2021).
While previous research has highlighted the importance of the variability introduced by PTFs in hydrological and land surface responses (e.g., Chirico et al., 2010;Mohajerani et al., 2021), there has been no comprehensive analysis on how those uncertainties could propagate from the water cycle to the carbon cycle via uncertainties in the predicted plant water stress, which limits plant water uptake and photosynthesis. Specifically, little attention has been given on how different shapes of the water retention curve, as well as differences in hydraulic conductivity, affect soil water flow and, subsequently, plant water availability and vegetation dynamics. Additionally, most land surface and ecosystem model applications, which are essential for climate change projections, neglect the role of topographic features as they consider one-dimensional vertical domains, thus potentially misrepresenting the importance of soil hydraulic properties (see also Fatichi et al. (2020)).
In this study, by using a state-of-the-art ecohydrological/terrestrial biosphere model with a dynamic vegetation component, we quantify how disagreements between PTF predictions lead to uncertainties in the dynamics of the water and carbon cycles at 79 locations worldwide, covering most terrestrial biomes and soil types. The analysis is first carried out for one-dimensional domains to mimic the common use of these models at global or large scales. Subsequently, we take small scale topography explicitly into account by simulating a distributed domain of 4 km 2 with a 40 × 40 m 2 resolution. The specific hypotheses of the study are: Numerical simulations are carried out with the ecohydrological and terrestrial biosphere model T&C (e.g., Fatichi et al., 2012Fatichi et al., , 2021Fatichi, Manzoni, et al., 2019;Mastrotheodoros et al., 2020;Paschalis et al., 2015Paschalis et al., , 2018Z. Zhang et al., 2021). T&C computes the water, carbon and energy balances at the land surface, and can be run as a distributed model on a regular grid geometry, taking into account small scale topography.
Specifically, in T&C, soil water transport is simulated using a quasi-3-D formulation of the Richards equation.
Vertical soil water fluxes are simulated by solving the 1-D Richards equation using a finite volume approach. Horizontal soil water fluxes are approximated using a kinematic wave approximation (Fatichi et al., 2012;Hopp et al., 2016). Surface water flow, including both hillslope and channel routing, is simulated using the kinematic wave approximation with a time adaptive finite volume discretization and a D8 downhill flow direction.
Evapotranspiration is simulated using the common resistor analogue scheme (Bonan, 2019). Resistances that are considered in T&C include atmospheric, leaf boundary, undercanopy, and soil resistances. Soil resistance affects soil evaporation and it is simulated using the formulation of Shahraeeni et al. (2012). Leaf stomatal resistance that affects plant transpiration is simulated using the Leuning-based formulation (Paschalis et al., 2017;Tuzet et al., 2003).
Soil water dynamics affect the carbon cycle through their impact on limiting photosynthesis during periods of low plant water availability. Photosynthetic rates, modeled according to Collatz et al. (1991Collatz et al. ( , 1992, are limited when the root zone integrated soil water potential drops below a plant-specific threshold. This threshold is a flexible parameter in T&C and can be tuned to represent any vegetation species. Photosynthesis suppression with decreasing water potentials follows a sigmoidal function (Fatichi et al., 2012). Prolonged duration of high vegetation stress can trigger leaf shedding in T&C, based on a plant specific stress duration threshold. When photosynthesis is limited, plant stomata close and stomatal conductance (i.e., the reciprocal of resistance) is reduced linearly with the reduction of leaf-scale photosynthesis. Vegetation structure affects hydrological fluxes via the dependence of interception, throughfall, and transpiration on leaf area.
The carbon assimilated through photosynthesis is either used for autotrophic respiration, which is temperature dependent, or used for plant growth by being allocated in the model's plant biomass pools (leaves, living sapwood, fine roots, carbohydrate reserves, fruits and flowers, and heartwood). Plant phenology is simulated using empirical formulae, where soil water stress, root zone temperature, solar radiation, and photoperiod are the factors leading to changes in the plant's phenological state.
In T&C both the BC and MvG soil hydraulic functions can be used. The modification of Webb (2000) is implemented for the MvG model to improve its predictions at low values of soil water content. In T&C the saturated hydraulic conductivity and the parameters describing the shape of the water retention curve can be derived by PTFs. In this study, seven different PTFs were used (Figure 1, Figure S2 in Supporting Information S1): three to estimate the parameters for the BC model (Cosby et al., 1984;K. Saxton et al., 1986;K. E. Saxton & Rawls, 2006) and four for the MvG model (Weynants et al., 2009;Wösten et al., 1999;Y. Zhang & Schaap, 2017;Zacharias & Wessolek, 2007). All PTFs, apart from that of Zacharias and Wessolek (2007), provide values for both the saturated hydraulic conductivity and the parameters for the shape of the water retention curve. For the case of Zacharias and Wessolek (2007) PTF, estimates for the saturated hydraulic conductivity were obtained by the PTF of K. E. Saxton and Rawls (2006). All PFTs depend on a subset of the following soil properties: percentage of clay, sand, silt, organic material, and bulk density.
Processes related to soil biogeochemistry and belowground carbon, nitrogen, phosphorus and potassium cycles are not simulated in this study, even though the model we used has this capability (Fatichi, Manzoni, et al., 2019).

Plot Scale Simulations
The first set of simulations performed here neglects topography and uses T&C at the plot scale, for example, a one-dimensional vertical model over a flat domain, similarly to global scale applications of Earth system models. The model was set up for 79 sites globally, covering a wide range of climates and biomes (Table S1 in Supporting Information S1). For these sites, the model has been previously calibrated and validated (Fatichi, Leuzinger, et al., 2016;Fatichi & Pappas, 2017;Moustakis et al., 2022;Paschalis et al., 2018). To run the model, we used hourly local weather data including air temperature, short and longwave radiation, humidity, cloudiness, and wind speed. The time span of data availability was between 2 and 32 years, with an average simulation length of 9.3 years.
For all 79 sites, two sets of simulations were carried out. For the first simulation set, which was used for the quantification of uncertainties due to PTF differences in ecosystem functioning, local soil properties (soil texture, organic content, and bulk density) were used, and the model was run using all seven PTFs. The second simulation set was used to compare the within soil-type uncertainties that are introduced by the choice of PTFs with the across soil-type uncertainties that are introduced by lack of knowledge of the exact soil physical properties. To this purpose, for each site 12 × 7 simulations were performed where 12 different soil types, representative of each USDA class, were used (Table 1), together with all 7 PTFs, resulting in a total of 6,636 simulations. In all soil textures, we added a 2% fraction of soil organic matter, a texture not included in the USDA classification. In this sense, our definition for each class is not identical to USDA but very similar (classification details are presented in Table 1). For the PTFs that did not include an organic matter fraction in the soil texture, this 2% fraction was considered as silt.

Distributed Simulations
For two of the 79 sites, Konza and Oak Ridge (ORNL), representative of a water-limited grassland and a temperate deciduous forest, respectively, the distributed version of T&C was used to quantify how spatially distributed water flow dynamics in a complex topography alters PTF-related uncertainties. Specifically, the model was run for both sites using an idealized 4 km 2 steep catchment topography ( Figure 2) with a 290 m elevation difference draining to a single outlet. Land use was considered identical for the entire domain (i.e., no spatial heterogeneity). All other properties (e.g., plant  (2007) (Ivanov et al., 2008;Remondi et al., 2019). The spatial discretization of the model was 40 m, leading to 2,500 computational cells.
For both sites, two soil types were used. One representative of a fast draining sandy soil (sandy soil: 40% sand, 5% clay, 2% organic, and a bulk density of 1.47 gcm −3 ) and one of a slow draining clay soil (10% sand, 30% clay, 2% organic, and a bulk density of 1.34 gcm −3 ). The soil depth was identical to the soil depth used for the point scale simulations for those two sites (Konza: 2,500 mm, ORNL: 1,500 mm). For the two sites and the two soil types, a 4-year long simulation was performed using all 7 PTFs. The soil type and soil depth were considered homogeneous across the entire catchment and an anisotropy ratio of 50 between the horizontal and vertical hydraulic conductivities was considered for all simulations (e.g., Assouline & Or, 2006). An additional simulation with an anisotropy ratio of 10 (not shown here) was also performed to assess the sensitivity of our simulations to this parameter, which turned out to be of low importance and thus not further discussed. A zero water flux boundary condition at the soil bottom was used in the spatially distributed simulations to avoid additional uncertainties related to resolving bedrock groundwater dynamics, by neglecting flow at the soil-bedrock interface. Accompanying the spatially distributed simulations, point scale simulations with identical meteorological forcing and soil properties (including the zero flux boundary condition) were performed for those two sites to facilitate comparison and allow a proper quantification of the sensitivity of accounting for topography and distributed water fluxes.

Uncertainty Analysis in Point Scale Simulations
For the first set of point scale simulations (i.e., those that used the observed local soil texture) and for the distributed simulations, we computed the uncertainties that PTF differences introduce in other variables as the average, across all time steps, of the standard deviation of modeled variables using different PTFs-that is, for a variable ℎ ( ) , where i is the PTF, h the time scale of aggregation, and t the time step, the uncertainties are computed as: where n T is the total number of time steps, n P the number of PTFs and ℎ ( ) To obtain nondimensional coefficients of variation (CV) at the annual scale, the average value of annual standard deviations was normalized with the mean annual value of the corresponding flux for all PTFs. The variables chosen to describe water fluxes were the evapotranspiration (ET), plant transpiration (T), abiotic evaporation (E), comprising both ground evaporation and evaporation from interception, total runoff (R), and leakage to deeper soil layers (Lk), defined as the water flux at the bottom boundary of the simulation domain. To facilitate comparison of water fluxes across sites, ℎ was normalized by the site's average precipitation at the same time scale. The variables chosen to describe carbon fluxes were the gross and net primary productivity (GPP and NPP), that is, the gross photosynthetic rate, and the net carbon used to produce plant tissues, that is, GPP minus autotrophic respiration. The structure of vegetation was described by the leaf area index (LAI). The plant water stress was quantified by the vegetation stress factor β, which is the water stress reduction factor applied to photosynthesis when soil water potential is below a threshold. β ranges between 0 and 1, with 1 indicating unstressed vegetation and 0 complete suppression of photosynthesis. β is a sigmoidal function of soil water potential for values of water potential below a plant specific threshold (Wu et al., 2018). For the distributed simulations, the variation of LAI and ET for every grid point were also normalized by their average value across all time steps and PTFs for any given grid point, to obtain nondimensional coefficients of variation.
For the second set of simulations at the point scale (i.e., those that used 12 different soil classes), an approach similar to Paschalis et al. (2017) was used to partition variability into its natural component (i.e., stochastic variability) and its uncertainty due to different soil textures and different PTFs. For example, for a variable Y i,j,k (t), where i corresponds to the PTF, j to the soil texture, k to the site, and t is the time step, the natural variability is quantified by: ( 2) where I, J, K, N are the number of PTFs, soil textures, sites and years accordingly. Similarly, the PTF uncertainty is quantified by And the texture uncertainty is: The variance partitioning was performed at the annual scale.
Uncertainty dependence on background climate was investigated using the wetness index of each site, defined as the average ratio of annual precipitation to annual potential evapotranspiration for all years in the meteorological record = .

The Role of Topography
To understand how consideration of topography alters uncertainties, we computed the uncertainties introduced by PTFs for all grid points within the distributed domain, for ET, LAI, root zone soil moisture, and β. The dependence of the uncertainties on the topography was computed as a function of each grid point's topographic index, defined as TI = log(A/tan(b)), where A [m 2 ] is the area draining to any given grid point and b [rad] the slope at the same grid point defined as the steepest downhill slope based on the D8 flow direction method (e.g., Wilson et al., 2007). The catchment response was quantified using the probability distribution of discharge, expressed in terms of its survival function (i.e., probability of exceedance) at its outlet at the daily time scale.

Plot Scale Simulations
The uncertainties introduced by the PTF disagreement in predicting soil hydraulic parameters are significant for all ecosystem variables. Considering annual scale statistics (Figure 3), uncertainty in vegetation structure is on average ( ) = 0.16 m 2 m −2 (Figure 3a) which is 8% of the LAI mean ( Figure S3 in Supporting Information S1). The uncertainties in GPP and NPP are ( ) = 51.4 gC m −2 y −1 and ( ) = 30.2 gC m −2 y −1 (Figure 3b), which are 7.6% and 9.8% of their mean values, respectively ( Figure S3 in Supporting Information S1).

of 18
The very similar magnitudes of LAI, GPP, and NPP uncertainties are expected as GPP directly depends on LAI, with photosynthetic rates scaling with LAI, and LAI in turn depends on GPP as it is controlled by the available carbon. As autotrophic respiration is less affected by plant water availability, the uncertainties of GPP and NPP are expected to be similar. GPP is affected by soil water availability (and thus PTF choice) in two ways: first, by short-term plant physiological responses that lead to stomatal closure and photosynthesis reduction, and second, by longer-term changes of LAI due to water availability.
Uncertainties introduced by PTFs in water fluxes at the annual scale were on average 4%, 3.4%, and 7.2% of the received precipitation, for evapotranspiration, runoff, and leakage, respectively ( Figure 3c). Even though uncertainties are a small fraction of the total water received as precipitation at each site, when water fluxes are compared with their annual averages, uncertainties were on average 5.8%, 49.8%, and 169.8% of their mean values for evapotranspiration, leakage, and runoff, respectively ( Figure S3 in Supporting Information S1). Runoff production is the most impacted variable in the experiments. Our results suggest that choosing among different PTFs can lead to large uncertainties with regard to runoff generation, without changing substantially the overall water budget. The reason for this behavior is that, in point scale experiments, simulated runoff is commonly a relatively small flux, as, by neglecting topography, we neglect flow convergence zones and thus the main locations for runoff generation. In addition, assuming a flat surface leads to infiltration reaching its maximum value, compared to steep hill slopes where infiltration rates may be lower.
The uncertainties of both hydrological and vegetation dynamics are dependent on the site's climate, albeit with a noticeable scatter (Figure 4). On average, uncertainties in evapotranspiration and plant productivity (NPP) increase in drier climates, the main reason being that, in arid/semiarid climates, vegetation is water stressed during part of the year, affecting both plant transpiration, which is the largest fraction of ET (Fatichi, Leuzinger, et al., 2016;Paschalis et al., 2018), and photosynthesis. In those climates, differences across PTFs lead to differences in the dynamics of soil water potential, thus affecting the strength and duration of plant water stress. In wetter climates, where water limitations occur rarely, uncertainties regarding productivity are small, despite considerable effects on soil water dynamics of the PTF choice. For climates where rainfall exceeds potential evapotranspiration (i.e., wetness index WI = P/PET > 1) uncertainties are mostly negligible, as plants do not experience significant water stress.
The uncertainties related to leakage (i.e., a proxy of groundwater recharge) are largest in intermediate climates (Figure 4). In arid/semiarid climates, irrespective of the PTF that describes the soil hydraulic conductivity and the shape of the characteristic curve, leakage is low due to persistent low soil moisture conditions. In mesic sites, where precipitation far exceeds potential evapotranspiration, uncertainties related to leakage are also small. In such sites, soil is commonly close to field capacity, and precipitation is mostly driving the amount of water recharge regardless of the soil type, with differences arising primarily due to PTF disagreement in the values of the saturated hydraulic conductivity. At intermediate sites, where soil moisture is fluctuating between dry and wet conditions, all uncertainties due to the saturated and unsaturated hydraulic conductivity as well as the shape of the characteristic curve can affect the simulated leakage. Runoff uncertainties are independent of climate, as they occur equally at all sites.
As many of the sites analyzed in this study experience a highly seasonal climate, uncertainty analysis was also performed at sub-seasonal scales, splitting seasons into dry and wet months, with dry months being defined as those when monthly potential evapotranspiration exceeds monthly rainfall and wet months the opposite. The uncertainty magnitudes and their dependence on climate are similar between dry and wet seasons ( Figure S12 in Supporting Information S1), with the only exception being plant productivity. Uncertainty in plant productivity during dry seasons is on average lower than in the wet seasons, primarily in semiarid sites, because plant activity is at minimal rates no matter the choice of PTF.

Soil Texture and PTF Choice
In Figure 5, we compare the total variability of both annual water and carbon fluxes due to the natural variability of climate (i.e., variability from 1 year to another due to random weather fluctuations) with uncertainties induced by soil texture and uncertainties due to the PTF choice. The interannual variability for NPP (and similarly GPP and LAI), ET (and similarly E and T) and Lk far exceeds the uncertainties due to either soil texture or the variability amongst PTFs (Figures 5a, 5b and 5d). This clearly illustrates the dominance of the model's forcing in a physically based model rather than its parameters in simulating these fluxes. What is of major importance, however, is that the uncertainties due to soil texture are comparable to, and in most cases lower than, the uncertainties across PTFs. This means that a complete lack of knowledge of the soil texture at a site introduces less uncertainty in the simulated water and carbon fluxes than the uncertainties introduced by the choice of a given PTF once the soil texture is perfectly known. The within-soil type PTF variability for both water and carbon fluxes is mostly soil-type independent, with the exception of silt and silt-clay-loam experiencing the larger uncertainties for evapotranspiration and leakage ( Figure S5 in Supporting Information S1) compared to the rest of the soil types used in this study.
Noticeably, the runoff uncertainties introduced by the stochastic variability of weather are similar in magnitude to the uncertainties introduced by soil texture and lower than PTF choice (Figure 5c). It should be highlighted, however, that a large part of this uncertainty relates to very low saturated hydraulic conductivity values predicted by the PTF of Weynants et al. (2009), leading to much higher values of runoff for this specific PTF. This disagreement is likely due to methodological differences in the quantification of saturated hydraulic conductivity. In the PTF of Weynants et al. (2009), the saturated hydraulic conductivity values used for PTF parameterization were computed using a parametric fit for conductivity values measured below saturation. This estimate eliminates the potential influence of macropore flow that would be present in all other PTFs, and for this reason, the estimates of conductivity of Weynants et al. (2009) are significantly lower (see also Figure 1a).

Distributed Simulations
Uncertainties introduced by PTFs in both hydrological and vegetation dynamics can be either amplified or dampened when the effects of topography through radiation and lateral water fluxes are taken into account. As the patterns between GPP and NPP are very similar to LAI, only results pertaining to LAI are presented here. However, these are also representative of both gross and net primary productivity.
In the temperate ORNL site, where vegetation is not significantly water stressed, uncertainties in vegetation structure and ET at the annual scale are always <10% compared to their mean values (Figures 6a and 6b). LAI uncertainties always decline with increasing TI for both high and low permeability soils and are dampened when compared to the point scale simulations. However, large uncertainties occur in the small fraction of the catchment that experiences frequent inundation. The reason for this is that, in T&C, photosynthesis is suppressed for partially submerged or submerged vegetation, and small changes in the time a grid cell remains inundated can lead to large uncertainties in vegetation productivity and, ultimately, LAI. Also, radiation fluxes in (partially) submerged vegetation and vegetation dynamics are treated in a simplified manner in T&C. For this reason, we refrain from further interpreting this result, as it is dependent on the spatial resolution used for solving surface water fluxes and on the implementation of water-logging stress, which has not been fully tested in T&C. Additionally, as T&C adopts the simplified kinematic wave approximation for surface water dynamics and not the full dynamic wave approximation of the shallow water equations, inundation extent, which is the largest source of uncertainty for very high TI values, is unlikely to be well captured. Uncertainties regarding ET in ORNL are amplified, compared to the point scale simulations, for the low conductivity clayey soil, with the largest amplifications occurring in the flow convergence zones (high TI) (Figure 6b). In contrast, uncertainties are dampened when a sandy soil was used, with the dampening effect being strongest at intermediate TI values, which correspond to hillslopes near, but not on, the river network where water flow converges. Within the river network, a large scatter was obtained, primarily due to transpiration suppression due to inundation. For the same reasons as for LAI, we refrain from further interpretation of the uncertainties in channel areas.
In the water-limited Konza site, the uncertainties introduced by PTFs in simulated LAI values are large (up to 30% of the mean LAI values). Compared to the point scale simulation, uncertainties are amplified when a clay-rich soil was used and dampened in the case of a sandy soil (Figure 6c). ET patterns follow a similar trend (Figure 6d). Clay soils always exhibit higher uncertainties in the distributed simulations than sandy soils, even though for the equivalent plot scale simulations the differences were minimal (Figure 6). Uncertainties in the clay-soil simulation were always higher (for the entire simulation domain) than the equivalent point scale results for the water-limited Konza case, while for the sandy soils, uncertainties were dampened upland and amplified close to the flow convergence zones. Within the Konza catchment, uncertainties for both LAI and ET increase at high TIs close to the river network, with an exception being LAI uncertainties for TI < 9 in sandy soils.
In the water-limited site, most of the topographic variability of LAI can be explained by the uncertainties of simulated plant water stress (Figures 7c and 7d). The pattern for ET follows an almost identical pattern as for LAI ( Figure S6 in Supporting Information S1), the reason being the dependence of transpiration, that constitutes the largest fraction of ET, on the amount of leaf area and thus NPP. Uncertainty in the stress factor β is not highly correlated with the uncertainties in soil moisture ( Figures S7 and S8 in Supporting Information S1). The reason is that β depends on soil water potential in T&C and thus the uncertainties of soil water content are modified by the heavily nonlinear relationship between water content and soil water potential. Uncertainties in water potential are typically larger ( Figure S2 in Supporting Information S1) at the drier end of the water retention curve, where plant water stress occurs. Overall, the soil hydraulic properties, the temporal structure of rainfall, and topography Conversely, in ORNL the impact of water stress on ET variability is mostly negligible, leading to lower total uncertainty values compared to the water-limited Konza site. In ORNL, both LAI and β uncertainties decline with the increase of TI for the sandy soil simulation, and are almost independent of TI for the clay soil simulation ( Figure S9 in Supporting Information S1). The pattern of ET cannot be explained by β (Figures S8 in Supporting Information S1), suggesting that most of the variability comes from differences in net radiation and surface soil evaporation which is influenced by surface soil moisture. Surface soil moisture is simultaneously impacted by soil hydraulic properties and by radiation availability, modulated by topography due to hill slope aspect and shading. The overall complex interactions between these processes shape the dependence of ET uncertainties on local topography.
One variable that is highly impacted by the differences across PTFs is river discharge (Figure 8), and particularly its extremes. Uncertainties for extreme events can be within an order of magnitude and occur in both sites and for both soil types. In the highly permeable sandy soils, discharge, both in terms of base flow and extremes, is lower than in clay soils. Using the PTF of Weynants et al. (2009) leads to the highest simulated values of discharge. The PTF of Wösten et al. (1999) also leads to high discharge extremes. For the simulations using the low permeability clay soil the PTFs of Wösten et al. (1999); Zacharias and Wessolek (2007); Weynants et al. (2009) lead to the highest discharge extremes. As expected, in all cases, high runoff production was associated with PTFs having the lowest hydraulic conductivity values (Figure S10 in Supporting Information S1; first boxplot), which remarks that when the conductivity values are low, there is a higher than average discharge simulated.
All the results presented in this section are generalizable across time scales, from hourly to yearly, as the same patterns occur at all scales. However, its worth noticing that uncertainties overall decrease with an increasing temporal scale of aggregation ( Figure S11 in Supporting Information S1).

Impact of Differences Among PTFs
While many studies have focused on differences across PTFs and their ability of accurately predicting soil hydraulic properties (e.g., Minasny et al., 1999), how uncertainties in soil hydraulic properties impact ecosystem dynamics at various spatial scales has received far less attention. Our results show that both hydrological and vegetation dynamics can be substantially affected by this source of uncertainty and that PTF choice might be as important, or generally more, than the choice of the soil type, as already found by Gutmann and Small (2007). Previous studies focused primarily on hydrological dynamics and our results are in agreement with their findings. Weihermüller et al. (2021) showed that differences across PTFs can lead to a coefficient of variation for ET of ∼5%, in line with our 5.8% estimate, despite the very different modeling approaches used. Interestingly, Weihermüller et al. (2021) report a large number of cases where some PTFs led to numerical instabilities in the HYDRUS model. This problem did not occur in the T&C simulations, partly because of the different numerical scheme (finite volumes vs. finite elements) and partly because of the correction of Webb (2000) we introduced in the MvG model to improve representation of the water retention curves in the driest regime. Our results were also in line with Edouard et al. (2018); Yang et al. (2018); Rieger et al. (2010); Mohajerani et al. (2021) showing large uncertainties regarding runoff production, even though we cannot directly compare the uncertainties due to differences in climates and topography with the case studies presented in previous research. However, a crucial result is that uncertainties in runoff generation are of the same order of magnitude in all climates and soil types analyzed ( Figure 4). As knowledge of precise hydraulic properties is not common in most parts of the world, this finding highlights the long-standing challenge in hydrology of making flood risk simulations based solely on a process-based approach, which entails irreducible uncertainties.
Beyond runoff generation, the hydrological flux mostly affected by uncertainty was water percolation to deeper soil layers (i.e., leakage), which is a proxy of groundwater recharge. Uncertainties regarding leakage were largest, as a fraction of precipitation, at locations of intermediate wetness. In such areas, pumped groundwater is commonly a large fraction of the total water use (Gleeson et al., 2012), highlighting the importance of the uncertainty propagation from PTF estimates to recharge dynamics. The total uncertainties could also be larger as PTFs are mostly trained on surface soil data (e.g., Rawls et al., 1991). To better tackle this problem, more data are needed regarding soil hydraulic properties of deeper layers to further parameterize PTFs.

Implications for Carbon Cycle and Terrestrial Biosphere Models
This study is the first to estimate how uncertainties introduced by PTF choice impact carbon cycle dynamics. Our results suggest that, on average, variability of soil hydraulic parameters can lead to ∼10% uncertainty in both carbon fluxes and vegetation structure (i.e., LAI). This uncertainty, even though large, is much smaller than the uncertainties reported between different terrestrial biosphere models (typically in the order of 10%-40% for GPP and NPP for different parts of the world; e.g., Keenan & Williams, 2018;Zheng et al., 2020;Cramer et al., 1999). Combined uncertainties related to differences in how they simulate processes from photosynthesis (e.g., Pappas et al., 2013;Rogers, 2014), to phenology (e.g., De Kauwe et al., 2017Richardson et al., 2013), carbon allocation (e.g., Fatichi, Pappas, et al., 2019;Franklin et al., 2012), and vegetation water stress (e.g., Wu et al., 2018;Paschalis et al., 2020) largely exceed the estimated uncertainties introduced by PTFs. However, an uncertainty for LAI and carbon fluxes in the order of ∼10% is still significant and comparable to changes expected with climate change (e.g., Shao et al., 2013). Overall, reducing epistemic uncertainties between current generation terrestrial biosphere models and the uncertainty introduced by their most sensitive parameters is still of high priority for better carbon cycle predictions (e.g., Lovenduski & Bonan, 2017;Pappas et al., 2013Pappas et al., , 2016. Even though PTF uncertainties introduced in LAI and carbon fluxes are less than other sources, their importance is still high, particularly in water-limited ecosystems. Such ecosystems are amongst the most vulnerable under a changing climate (e.g., Huang et al., 2017), and projected increases in atmospheric, meteorological and hydrologic drought can intensify plant water stress. Changes in plant water stress can alter both the local water and carbon cycles, with pronounced impacts on global scale carbon dynamics (e.g., Humphrey et al., 2018) and also on the livelihood of people, as in many semiarid areas globally, subsistence is to a large degree dependent on rainfed agriculture, particularly in the developing world. It is fundamental to remark that a proper quantification of plant water stress cannot ignore a robust definition of soil hydraulic properties. While a large number of recent studies have focused on the importance of representing plant hydraulic and physiological properties in simulating water stress in terrestrial biosphere models (e.g., Feng, 2020;Kennedy et al., 2019;Medlyn et al., 2016;Xu et al., 2016), much less attention has been dedicated to soil hydraulic properties, which might be equally important in defining plant-water stress.
PTFs-induced uncertainties in plant water stress were able to explain, to a large extent, the uncertainties in both water and carbon dynamics. However, uncertainties in plant water stress were not highly correlated with uncertainties in soil moisture ( Figure S8 in Supporting Information S1). The reason is that, in the model we used, plant water stress depends on root integrated soil water potential (which is the quantity felt by plants) and not soil moisture. Many alternative models use soil moisture for parameterizing plant water stress (e.g., Harper et al., 2021;Paschalis et al., 2020;Wu et al., 2018), and thus the uncertainties of both the water and carbon dynamics, particularly in drier areas, are expected to be less PTF-dependent and more dependent on how soil moisture is used to parameterize water stress. This highlights the importance of simultaneously monitoring both soil water content and vegetation stress dynamics, in order to reduce uncertainties and, at the same time, improve and best tune terrestrial biosphere models. New generation remote sensing technologies can now provide to some degree such information at the global scale (e.g., Chan et al., 2016;Konings et al., 2021;Y. Sun et al., 2015).

The Importance of Accounting for Topographic Effects
Small scale topography was found to both amplify and dampen (based on soil type and climate) uncertainties introduced by the choice of the PTF. Particularly, in dry conditions, uncertainties within the spatially distributed simulations were up to threefold larger than plot scale simulations, for both water and carbon dynamics. This effect is largely related to lateral water redistribution in the surface (runoff and runon processes) and subsurface, which is considerably affected by soil hydraulic properties. In a distributed domain, such effects are accumulated from upland areas toward down slope converging areas, considerably increasing the importance of soil hydraulic properties in the simulation of plant water stress and thus carbon and water fluxes. To a second order, topographic variability also implies a different amount of received net radiation at the surface, which interacts with ground evaporation. Such a process is also mediated by soil hydraulic properties. Currently, all global scale applications of ecosystem models, either offline or coupled with atmospheric models, neglect subgrid topography and lateral water fluxes. This potentially leads to an underestimation of the importance of properly representing soil hydraulic properties Tafasca et al., 2020). This might have implications for a number of dynamics, such as land surface coupling (Gevaert et al., 2018), runoff predictions (e.g., Ukkola et al., 2018;X. Zhang et al., 2014) and predictions of plant water stress and its impacts on the carbon cycle (e.g., Anderegg et al., 2013;De Kauwe et al., 2020). Introduction of new generation hyper-resolution models can potentially resolve these problems, albeit with large computational costs involved (e.g., Fan et al., 2019;Forrester & Maxwell, 2020;Mastrotheodoros et al., 2020;O'Neill et al., 2020).
One additional source of uncertainty related to soil hydraulic properties, not investigated here, is the role of soil structure (Bonetti et al., 2021;Fatichi et al., 2020). All PTFs used in this study neglect the effect of soil structure that can significantly alter the soil hydraulic parameters at water potentials close to saturation. In fact, the uncertainties expected from soil structure alone have been found to have a similar order of magnitude as the one reported in this study, particularly for hydrological dynamics . Ultimately, uncertainties introduced by soil structure are further increased by complex soil-vegetation combinations (Bonetti et al., 2021), and can significantly impact hydrological and land surface dynamics. Based on results obtained here, it is likely that the role of soil structure in modifying ecohydrological variables will be magnified by considering a distributed topography, as already hinted by Fatichi et al. (2020).

Study Limitations
The results of the present study were derived using a single ecohydrological model. One source of uncertainty that was not considered here is the role of the model structure itself. The different estimates of soil hydraulic properties across PTFs should affect soil water transport similarly amongst most ecohydrological and terrestrial biosphere models, all other things being equal. That is because most models solve the Richards' equation (albeit exceptions that use simplified hydrological dynamics exist; see Wu et al. (2018) for a review of how different models solve soil water dynamics), commonly using similar numerical procedures. However, the representation of how soil moisture impacts plant water stress and thus plant transpiration and productivity is highly variable across models (e.g., Egea et al., 2011;Paschalis et al., 2020;Wu et al., 2018;Zhou et al., 2013). These model-specific differences will impact all ecohydrological dynamics. Similar large differences across models can be expected for the parametrization of the dependence of soil evaporation on surface soil moisture (e.g., Lehmann et al., 2018Lehmann et al., , 2008 and also how plants access water in the root zone (e.g., Harper et al., 2021). To compute the magnitude of uncertainties related to the choice of the ecosystem model, a multimodel intercomparison project would be necessary, or even better, an analysis where single model components are replaced with alternative solutions within the same modeling framework.
In this study, we analyzed how uncertainties of soil hydraulic properties propagate to hydrological and ecosystem dynamics. However, disentangling the exact mechanism of how those uncertainties propagate across spatial and temporal scales is a major challenge. Many of the variables we analyzed are highly correlated ( Figure S13 in Supporting Information S1) and causally dependent. For example, uncertainties regarding plant water stress will, at short time scales, impact plant transpiration and GPP/NPP through stomata closure. At longer time scales, LAI is reduced if GPP decreases, and it can subsequently impact both water and carbon fluxes. To disentangle the exact mechanisms, multiscale causal inference techniques would be needed (e.g., Detto et al., 2012;Runge et al., 2019).
Finally, we used for all sites and soil texture combinations all PTFs, without filtering out potentially unrealistic results, as the scope of the study was to compute the uncertainties originated by the indiscriminate use of PTFs, beyond specific areas and limits of derivation as often done in large-scale models. Developing a filtering scheme that could identify the applicability range of each PTF would likely reduce their uncertainties. In this sense, our estimates are likely conservative and quantify the maximum values of the PTF related uncertainties. Approaches that can successfully reduce uncertainties related to PTFs have been previously applied (e.g., McNeall et al., 2016;Verhoef & Egea, 2014), and could be incorporated in global scale modeling applications where currently PTFs are used without any filtering.

Conclusions
In this work, we quantified how disagreement between predicted hydraulic properties by 7 common pedotransfer functions leads to uncertainties in the simulated carbon and water fluxes at 79 sites across the world. We found that PTF disagreement can significantly impact both the water and carbon cycles, irrespective of soil type. Runoff and leakage to deep soil layers were the two hydrological fluxes most impacted by this uncertainty, with the uncertainty of annual fluxes being on average 168% and 50% compared to their mean values. Runoff uncertainties were independent of climate, while leakage uncertainties were strongest at climates of intermediate wetness (i.e., a wetness index of around 1). Despite those large uncertainties, because of the relatively small magnitude of those fluxes compared to annual rainfall, the uncertainties they introduce on the partitioning of precipitation to evapotranspiration and the other hydrological budget components is within 10%. Uncertainties in leaf area index and carbon fluxes (GPP and NPP), computed as the annual coefficient of variation, were smaller (8%, 7.6%, and 9.8%) than water fluxes and to a large extent explained by the intensity and duration of plant water stress. For both hydrological and vegetation dynamics, the uncertainties due to PTF disagreements were larger than the uncertainties due to differences in soil types, highlighting the importance of the choice of a PTF and the need for new, more robust and generally applicable PTFs. Finally, the results showed that small scale topography can amplify the importance of the uncertainties of soil hydraulic parameterization by threefold, especially for low permeability, clay-rich soils in water-limited ecosystems, highlighting that new generation hyper-resolution and distributed climate and terrestrial biosphere models might be more affected by soil hydraulic parameters than the current ones.

Data Availability Statement
The model used for the simulations can be found in the following citable repository, https://doi.org/10.24433/ CO.0905087.v2.