Stochastic analysis of plant available water estimates and soil water balance components simulated by a hydrological model

The uncertainty in soil hydraulic parameters is often not taken into account in process‐based hydrological modeling. Performing runs with 104 stochastic parameter realizations, we evaluated the propagation of uncertainty in the Van Genuchten–Mualem (VGM) parameters into estimates of the threshold values of soil water content used to calculate the total and readily available water, and on the long‐term (30 years) simulations of evaporation, transpiration, bottom flux, and runoff by the SWAP hydrological model. The simulated scenarios included weather data from a location in southeast Brazil and seven soils from the same region cropped with maize, comprising a wide range of texture classes. The results showed that uncertainties in VGM parameters affect the estimates of total and readily available water. Water balance components obtained by a deterministic simulation with average VGM parameters did not always agree with the average or median of stochastic simulations, and stochastic simulations including parameter uncertainties should be preferred. Variations in yearly rainfall characteristics were more important for bottom flux and evaporation, while transpiration and runoff were more strongly influenced by the variations in soil hydraulic properties.


INTRODUCTION
Hydrological modeling is a tool for predicting the effect of soil water balance components on ecosystem services (Decsi et al., 2022), agricultural crop growth (Siad et al., 2019), groundwater recharge and pollution (Amin et al., 2017;Nolte et al., 2021), erosion, and soil loss (Shojaei et al., 2020), among other practical applications.Hydrological models based on the Richards equation require the quantitative description of fundamental soil hydraulic properties, namely, water retention and hydraulic conductivity functions.In a process-based approach, atmospheric conditions together with properties of the canopy-atmosphere boundary layer are used to determine the potential transpiration rate, whereas soil hydraulic properties and soil-root boundary conditions control the potential root water uptake (RWU) rate (Couvreur et al., 2012;De Jong Van Lier et al., 2013).If the potential RWU rate exceeds the potential transpiration rate, the scenario is sink-limited, that is, the actual transpiration equals the potential rate.On the other hand, if the potential RWU rate is lower than the potential transpiration rate, crop Vadose Zone Journal transpiration will be source-limited, determined by soil water status and root system properties.
In this context, the estimation of water availability is an essential part of soil water balance modeling, as it determines the fraction of soil water available for crop transpiration.This estimation can be performed using the threshold values of soil water content used to calculate plant available water, field capacity (FC), wilting point (WP), and limiting point (LP), commonly linked to cascade reservoir or "bucket" models (Švob et al., 2022).Such models mimic the involved processes without addressing the physical mechanisms and do not explicitly include the soil hydraulic properties (De Melo et al., 2023;Jarvis et al., 2022).
Several approaches are also available to predict the actual RWU and transpiration rates as a function of the soil water status.One of these, the RWU function by Feddes et al. (1978), is included in the vadose zone hydrological models Hydrus (Šimůnek et al., 2016) and SWAP (Kroes et al., 2017).According to the Feddes function, actual transpiration is determined based on empirical threshold pressure heads defining the onset of transpiration-limiting conditions.Respective values for several crops are compiled in Taylor and Ashcroft (1972).
Besides the Feddes function, the most recent versions of the SWAP model (versions 4.X) also include the process-based RWU function MFlux (De Jong Van Lier et al., 2013).This function uses an axisymmetric single-root scenario, which is upscaled per soil layer and allows predicting RWU and transpiration based on soil hydraulic properties, root system geometry, and internal plant resistances (Casaroli et al., 2010;De Melo & De Jong van Lier, 2021;Durigon et al., 2012;Pinheiro et al., 2017).
The equations for water retention and hydraulic conductivity functions are presented as a set of equations with soil-specific parameters (Madi et al., 2018).Parameters can be estimated using measured data of soil water retention (θh) or hydraulic conductivity (K-h or K-θ).If enough data are available, pedotransfer functions can be developed to establish correlations between soil hydraulic parameters and other soil properties, most commonly based on information about particle size distribution, organic matter content, and bulk density (Minasny & Hartemink, 2011;Schaap et al., 2023).
The uncertainty of the parameter estimates may be due to measurement errors, intrinsic variability between samples, or inadequacy of the mathematical equation to mimic the physical phenomenon.This parameter uncertainty and its propagation into simulation results is often not taken into account in process-based hydrological modeling (e.g., Pinheiro & de Jong van Lier, 2021;Wesseling et al., 2020).
The propagation of parameter uncertainty into simulation results cannot be analyzed by a one-at-a-time sensitivity analysis, as parameters are correlated.An adequate approach to assess hydrological model output uncertainty due to the

Core Ideas
• Hydraulic parameter uncertainty and correlations can be used in stochastic simulations with a hydrological model.uncertainty in soil hydraulic parameters is by stochastic methods.Specific fitting software like RETC (Van Genuchten et al., 1991) reports the parameter uncertainty in terms of a standard error.Furthermore, correlations between the parameters are quantified in a correlation matrix.This information can be used in a Monte Carlo simulation approach (Caflisch, 1998;Ross & Marshak, 1991).Whereas a deterministic simulation with the mean parameter values yields one single model output, a Monte Carlo (stochastic) approach results in a large number of output values, which are then analyzed by statistical methods to obtain, for example, means, medians, standard deviations, or uncertainties (e.g., A. K. B. Dos Santos et al., 2022;Pinheiro & de Jong van Lier, 2021).
In this study, we aimed to evaluate the uncertainty of flux-based threshold values for plant water availability and the uncertainty of the soil water balance components simulated by a hydrological model using a stochastic modeling approach.Specific objectives were to (1) conduct a stochastic analysis to evaluate the propagation of uncertainty of soil hydraulic properties into the threshold values of plant available water, FC, LP, WP; (2) investigate the impact of the uncertainty of these parameters on long-term (30 years) hydrological simulations of soil water balance components for seven southeast-Brazilian soils covering a wide range of texture classes; and (3) compare the stochastic to the deterministic estimates of plant available water and soil water balance components in some simulated scenarios.

Study area and soil sampling
Seven soils under agricultural (arable) use, labeled A-G, were sampled in the São Paulo state in southeast Brazil.At each location, undisturbed and disturbed samples were taken at about 10cm depth, representing the 0-20 cm layer, and at about 30-cm depth, representing 20-40 cm.Undisturbed samples were taken in large sample rings (8-cm diameter and 5-cm high, volume 250 cm 3 approx.),six replicas, as well as in small rings (5-cm diameter and 3-cm high, volume 60 cm 3 approx.),five replicas.Disturbed samples for particle size determination were taken at the same depths.The sampled soils cover a wide range of texture classes and three important USDA soil orders (Table 1).

Soil hydraulic properties
For each soil layer, three of the large undisturbed soil samples were saturated with deaired 0.005 M CaSO 4 solution and used to determine the saturated hydraulic conductivity (K s ) using a falling head permeameter (Ksat equipment by Meter Inc.).Measured K s values were corrected for temperature to 20˚C.The other three large samples were saturated with water and submitted to standard HyProp (Meter Inc.) analysis (Peters & Durner, 2008), in which two tensiometers were vertically aligned along the sample axis and placed at 1.25 and 3.75 cm from the bottom of the sample.The sample and tensiometers were placed on a balance and weights and tensions were automatically logged every 15 min, yielding a series of h-θ and h-K data pairs in the range of h between approximately −10 and −800 cm.The five smaller samples were saturated and subsequently submitted to pressure heads of −10, −20, −40, −60, −100, −330, and −600 cm on a tension table (EcoTech pF laboratory station) and pressure heads of −1000 and −3000 cm in a porous plate pressure chamber and weighed after equilibration (varying from a few days to a week depending on the tension).Disturbed samples were analyzed in a WP4 dewpoint tensiometer to yield water contents corresponding to pressure heads in the range between −5000 and −50,000 cm.Soil hydraulic properties were expressed as parameters of the Van Genuchten (1980) equations with the Mualem (1976) parametric restriction: and where S e is the effective saturation given by and h (cm) is the pressure head, θ (cm 3 cm −3 ) is the water content, θ s (cm 3 cm −3 ) is the saturated water content, and θ r (cm 3 cm −3 ) is the residual water content; α (cm −1 ), n, and l are empiric parameters; and K s (cm day −1 ) is the saturated hydraulic conductivity.The parameters of these equations will be referred to as the Van Genuchten-Mualem (VGM) parameters.
To ensure normality of the distribution of fitted parameters, K s , α, and n, improve the quality of the subsequently generated set of stochastic VGM parameter realizations, the transformations proposed by Carsel and Parrish (1988) were applied, that is, Inserting Equation (4) in Equations (1 and 2) yields, and Means, standard deviations, and the correlation matrix for parameters θ s , θ r , τ, ν, l, and κ were obtained by fitting Equations (5 and 6) to the measured data, that is, values of K s , h-K data pairs obtained from the HyProp equipment, values of total porosity (corresponding to θ s at h = 0), and h-θ data pairs obtained from the HyProp, the tension table, the pressure chamber, and the WP4 dewpoint tensiometer.Fitting was performed using specific software developed by Inforsato et al. (2022), reporting τ, ν, and κ, as well as l, and θ 1 and θ 2 (corresponding to h 1 and h 2 , taken in this study as h 1 = −1000 cm and h 2 = −10 cm), related to θ r and θ s by the equations as follows: ) This procedure considers θ r and θ s as mere fitting parameters, and negative values of θ r may occur without numerical consequences in subsequent simulations.

Generating stochastic realizations
The means and standard deviations of the transformed parameters along with their correlation matrix were used to generate 10 4 stochastic parameter realizations for each layer of the seven soils using a Cholesky decomposition performed by the software StochHyProp (De Jong van Lier, 2023), without tail exclusion.The resulting transformed parameters τ, ν, κ, θ 1 , and θ 2 were transformed back to their original VGM values (α, n, K s , θ r , and θ s ) using Equations (4, 7, and 8).Attempts to make an individual analysis of parameter values are generally ineffective as hydraulic parameters act together to result in a specific curve shape and hydrological behavior.The overall effect of parameters and their statistical properties is better evaluated concerning crop water availability or water balance components, as specified in the following items.

Evaluating threshold values for plant water availability
Threshold values for plant water availability, that is, FC, LP, and WP, were determined by flux-based methods considering a rooting depth of 60 cm.The calculations were performed for each of the 10 4 stochastic parameter realizations of each of the seven soils (A-G).The soil hydraulic properties below 20-cm depth were considered to correspond to the second measured layer (Table 1).
FC was determined using the SWAP hydrological model (Kroes et al., 2017) in a scenario similar to the one used by De Jong Van Lier (2017): a 120-cm deep soil profile was discretized into 24 layers of 5 cm.The initial hydraulic conditions were defined as corresponding to a saturated soil profile (pressure head h = −0.1 cm over the entire profile).Evaporation and precipitation were set to zero, and the lower boundary condition (at 120 cm depth) was set to the "free drainage of soil profile" option, that is, gravitational flow.Simulations were run for 90-day periods with an output every simulated hour.FC was considered to correspond to the hydraulic conditions at the first occurrence of a bottom flux (BF) <1 mm day −1 at the depth of 60 cm and expressed as the corresponding average θ or h between the surface and the 60-cm depth.
The LP, defining the onset of the falling rate phase of crop transpiration (Feddes & Raats, 2004;Metselaar & De Jong Van Lier, 2007), was calculated using the equation by Pinheiro et al. (2017), yielding the corresponding matric flux potential M LP (cm 2 day −1 ) given as follows: T A B L E 2 Mean Van Genuchten-Mualem (VGM) parameters (Equations 1 and 2) for both layers of the seven evaluated soils.

Soil
where K is the hydraulic conductivity as a function of the pressure head h, h ref = −1.5•10 4 cm, h LP (cm) is the pressure head at the LP, T p (cm day −1 ) is the potential transpiration taken as 0.4 cm day −1 , k is the number of layers (in our scenarios: k = 2), L i (cm) is the thickness of layer i (here: 20 and 40 cm, respectively), and Ρ i (cm cm −3 ) is the root length density per layer.We used Ρ 1 = 1.0 cm cm −3 and Ρ 2 = 0.1 cm cm −3 , which are common values for arable crops (see, e.g., the compilation of data by De Willigen & Van Noordwijk, 1987).Using these values, M LP = 2.817•10 −2 cm 2 day −1 .The corresponding value of pressure head h LP for each stochastic realization was found numerically using the first two terms of Equation ( 9), and the water content θ LP was calculated using Equation (1).The WP, corresponding to a low soil water content at which plants do not recover after 24 h in a zero-transpiration environment, was estimated using the matric-flux-based proposal by De Melo et al. (2023), using a multiplicative factor f WP : In this study, we used the value of f WP = 0.01 proposed by De Melo et al. (2023).It corresponds to M WP = 2.817•10 −4 cm 2 day −1 , and subsequently allowed us to calculate h WP and θ WP .
These threshold values for plant available water (FC, LP, and WP) expressed as soil water content (θ) were used to calculate the total available water, TAW (TAW = θ FC −θ WP ) and the readily available water, RAW (RAW = θ FC −θ LP ) contents.To express TAW and RAW in units of water depth (cm) for the rooted soil profile, they were multiplied by the rooting depth of 60 cm.

Crop simulation scenarios
Simulations were performed with the hydrological model SWAP, v. 4.2.32, an updated version of the fully documented v. 4.0.1 (Kroes et al., 2017).SWAP is a 1-D model that numerically solves a discretized version of the Richards equation with a root water extraction sink term (Van Dam & Feddes, 2000), describing vertical water flow in the soilplant-atmosphere environment.To do so, the model makes use of the VGM soil hydraulic parameterization.A rainfed maize crop was simulated using a simple crop file "MaizeS.crp"(Appendix A), including information about the rooting depth (z R ) as a function of the development stage reaching a maximum value of 100 cm, with a root length density of 1.0 cm cm −3 at the surface, linearly decreasing to zero at z R .The crop cycle was defined at a fixed length of 111 days, starting on December 1 of each simulated year and ending on March 21 of the subsequent year.
Hydrological simulations were performed for 30 years starting on July 1, 1992, and ending on June 30, 2022, encompassing 29 growing seasons starting on July 1 and ending on June 30 of the subsequent year.Weather data from the meteorological station of the University of São Paulo in Piracicaba, Brazil (22.7030˚S; 47.6235˚W, 526 m a.s.l), were used in the simulations (Table 1).Piracicaba is characterized by a winterdry tropical climate, Koeppen Aw.The annual average rainfall T A B L E 3 Average profile water content (0-60 cm) at field capacity (θ FC ), at the limiting point (θ LP ), at the wilting point (θ WP ), readily available water (RAW), and total available water (TAW) estimated deterministically with mean Van Genuchten-Mualem (VGM) parameters (VGM mean ) and stochastically as the mean estimate for 10 4 stochastic parameter realizations (mean r ) with the respective coefficient of variation (cv r = σ r /mean r ), for each of the evaluated soils. is 1300 mm, the average T min and T max are 15.4 and 28.1˚C, respectively, and the average radiation is 14.5 MJ m 2 day −1 .The monthly distribution of rainfall and temperature is shown in Figure 1.

Soil
To simulate transpiration and growth reduction due to water stress, the simulations were performed using the processbased RWU function MFlux (De Jong Van Lier et al., 2013), (SwDrought = 2 in the SWAP crop file-Appendix A).MFlux is based on the simulation of soil water flow toward individual roots.The corresponding axisymmetric water flow toward the roots is numerically solved using the matric flux potential, M (the integral of soil K[h] at a predefined h interval) as the driving variable (De Jong Van Lier et al., 2008, 2013).The concept also includes radial and axial hydraulic resistances within the plant system.
The MFlux function was used with root and xylem radii of 0.5 and 0.2 mm, respectively, a hydraulic conductance between leaf and root xylem L l = 1.03•10 −4 day −1 , and a radial hydraulic conductivity of the root tissue K root = 3.5•10 −5 cm day −1 .Other parameters from the MFlux function were set to their standard values, that is, root coefficient a = 0.53, root efficiency = 1.0, and the limiting leaf pressure head h l = −15,000 cm.

Generating simulated soil water balance components
The propagation of uncertainties in the soil hydraulic parameters and its effect on simulated soil water balance components were evaluated using the generated stochastic realizations of the VGM parameters (Section 2.3).For each soil (A-G), the SWAP model was run 10 4 times, once for each stochastic realization, yielding 10 4 realizations of the yearly soil water balance components bottom flux (BF), runoff (R), evaporation (E), and transpiration (T), for the period between July 1, 1992, and June 30, 2022.

Soil hydraulic properties
A summary of the parameter estimates is shown in Table 2, containing the backtransformed mean VGM parameters (Equations 1 and 2) for the seven soils.The full set of results from the fitting procedure is presented in Appendix A for both layers of the seven analyzed soils, including the Example of data points obtained by different methods and fitted curves (blue lines) for water content (Equation 1) and hydraulic conductivity (Equation 2) as a function of pressure head for one of the analyzed layers (soil F, 20-40 cm).
transformed parameter estimates with their standard deviations and the correlation matrix used to generate the stochastic parameter realizations.The l parameter has an average coefficient of variation of the order of 30%, reflecting the uncertainty in the measurements of unsaturated hydraulic conductivity.The coefficient of variation of the other parameters was around 3% on average, and an order of magnitude higher (around 10%) for ν = ln(n − 1).A very high negative correlation (<−0.95) was observed between parameters l and ν = ln(n − 1) in most of the analyzed soil layers.Parameters l and n are used in Equation ( 2) to calculate hydraulic conductivity K, and the observed correlation means that high values of l imply values of n close to 1.A high negative correlation was also obtained between l and τ = ln(α), and a high (positive or negative) correlation occurred between l and κ = ln(K s ) in some cases.
It should be noted that all parameter values are based on data obtained from the measurement techniques described in Section 2.2.The very high values of K s for the surface layer of soils B and D go together with values of n very close to 1, causing a very rapid decline of K with a decreasing pressure head close to saturation.In the case of the surface layer of soil  D, for example, the fitted value of K s = 5777 cm day −1 reduces to less than 1% of its value (K = 31 cm day −1 ) at h = −1 cm.This is unrealistic and shows that K s is a mere fitting parameter.On the other hand, when using these parameter values in a vadose zone hydrological simulation, the very high K s will only affect the results in fully saturated conditions.Figure 2 shows an example of data points obtained by the different methods and the fitted θ-h and K-h curves for one of the analyzed soil layers (soil F, 20-40 cm).In the example of Figure 2, but also in the other examined layers, differences between methods and variability due to replication can be observed.Variability of this order of magnitude may be caused by hysteresis, equipment calibration, or other sources (Gubiani et al., 2013;Guillaume et al., 2023;Schelle et al., 2013).It will determine the standard deviation of fitted parameters and the uncertainty of modeling results.Without a stochastic analysis, the impact of uncertainty on the simulated results would not be quantifiable.

Threshold values for plant water availability
For the seven soils, Table 3 shows θ FC , θ LP , θ WP , RAW, and TAW according to a deterministic estimation procedure using only the mean VGM parameters and the mean estimate for the 10 4 stochastic parameter realizations with the coefficient of variation.There are no important differences between the deterministic and the stochastic mean values, except for the case of soil D, where differences show up, especially in RAW and TAW estimates.The surface layer of soil D has a high uncertainty in α and K s , with a coefficient of variation for τ = ln(α) and κ = ln(K s ) both in the order of 10%, which in turn affects the coefficients of variation in this soil as well as the estimates of the threshold water contents.
Coefficients of variation for RAW and TAW are higher, for most of the soils, than the coefficients of variation for the three specific water contents (θ FC , θ LP , and θ WP ).As the calculation of RAW and TAW includes the difference between θ FC and θ LP or θ WP , the coefficient of variation will likely increase in RAW and TAW values.Figure 3 illustrates this with an example of θ LP versus θ FC estimated for the 10 4 stochastic realizations of soil A. Although the relative differences in the threshold water contents are small, the obtained data points cover a wide range of RAW, increasing threefold from about 1 to 3 cm, considering the 0-60 cm profile.
Figure 4 shows the cumulative frequency of RAW and TAW in soils D and G.In this example, the variability of values is higher in soil D, especially for RAW.Again, this is most likely caused by the higher uncertainty in α and K s in the surface layer of this soil, where both parameters affect the conversion of the limiting matric flux potential M LP (Equation 9) to the soil water content θ LP .
Given the stochastic realizations for each soil, no correlation was found between the threshold water contents θ FC , θ LP , and θ WP , on the one hand, and available water contents TAW and RAW, on the other hand.This indicates that the cumulative frequency of TAW and RAW, as shown in Figure 4, is not determined by one of the threshold water contents alone, but rather by both θ FC and θ WP in the case of TAW, and θ FC and θ LP in the case of RAW.This agrees with the fact that a positive correlation coefficient r (r ≈ 0.4-0.7)was found between θ FC and θ LP , and somewhat lower (r ≈ 0.2-0.5) between θ FC and θ WP for the stochastic realizations.Therefore, an increase in θ FC is likely to be accompanied by an increase in θ LP and θ WP , reducing the effect on TAW and RAW.

Simulated soil water balance components
Considering the 30 simulated years (1992-2022), the mean potential transpiration of the maize crop was 49.2 cm over the 111 days growing season, corresponding to an average of 0.443 cm day −1 , with a standard deviation of 8.6 cm over the years.The mean potential evaporation was 96.0 cm year −1 (equivalent to a daily average of 0.263 cm day −1 ), with a standard deviation of 7.1 cm.
The 30 evaluated years were organized in terms of the yearly (July 1-June 30) rainfall and the maize growing season (December 1-March 21) rainfall.About half of the yearly rainfall was normally observed during the maize growing season, from December to March (Figure 5).Considering both the yearly and the growing season rainfall, 2012-2013 showed to be an average year for both; the driest year was 2013-2014, while the wettest year was 1998-1999.These 3 years F I G U R E 6 Transpiration and bottom flux in soil A simulated for the 30 analyzed years in Piracicaba, Brazil, deterministically (solid lines) using the VGM mean parameter set and stochastically (dashed lines) using the mean value of the 10 4 parameter realizations.

F I G U R E 7
Cumulative frequency of simulated soil water balance components bottom flux, evaporation, transpiration, and runoff resulting from 10 4 stochastic realizations of hydraulic parameters of soil A and soil G for a dry (2013-2014), wet (1998-1999), and average (2012-2013) year.Dotted lines labeled "D" mark the deterministically obtained value.were labeled "average," "wet," and "dry," respectively, and analyzed with more detail in some examples in the following.
The average simulated values of the four soil water balance components, BF, E, T, and R, considering the 30 years, seven soils, and 10 4 stochastic realizations, resulted in the values shown in Table 4. Values of R are relatively insignificant for most soils and even zero for others and D).Values of E are very similar among the soils, but T showed greater differences and with higher standard deviations among years (σ y ).Consequently, BF tends to be lower in soils with higher T and vice versa.The standard deviations among the years are similar among the soils, with a coefficient of variation of about one-third of the average for BF and one-tenth for E and T. Comparing to the standard deviation among the soils Vadose Zone Journal (σ s ), σ y is three to five times greater for BF and E, showing the yearly rainfall characteristics to be more important for these components than the soil properties.In contrast, in the case of T and R, σ s > σ y and soil properties show to be more important.The median values (not shown) were very close to the average values, with differences <0.1% in all cases.
Figure 6 shows, for soil A, the transpiration and BF simulated along the 30 analyzed years, comparing the mean value from 10 4 simulations and the value from deterministic simulation based on the mean hydraulic parameter values (VGM mean , Appendix B, Table A1).In absolute values, the stochastically determined mean BF is always greater, as are the stochastic T values.Furthermore, higher variability of BF is observed compared to T. BF responds directly to rainfall with a smaller interference of soil hydraulic properties and crop water use, whereas T is the result of the strong interaction between rainfall, atmospheric demand, and the crop response to excess or deficit soil water conditions.
To analyze the results for BF and T in more detail, Table 5 contains the average and standard deviation of the 10 4 stochastic realizations for each of the seven soils simulated for the wet, average, and dry year, together with the deterministically obtained values using the VGM mean parameter set.Important differences between stochastic and deterministic simulation results can be identified.Stochastic BF is higher than the corresponding deterministic value in the wet year.In the average year, results are more similar, and in the dry year, stochastic BF is lower than deterministic BF.There is no straightforward explanation for these results, as they have their origin in complex parameter interactions in hydrological simulations.It is more important to observe that the parameter uncertainty, when considered using a stochastic approach, may lead to important differences in the simulated values of water balance components.Similarly, for the T component, differences between stochastic and deterministic values are observed, not so much in the wet year when T tends to its potential value with less influence of soil properties but pronounced in the dry year in some of the soils (C and E).
As an illustration of the results obtained from the stochastic analysis, Figure 7 shows cumulative frequency curves of BF, E, T, and R for the dry, average, and wet years in two of the soils, soil A representing soils with an increase in T from the wet to average to dry year and soil G representing those with the highest value of T in the average year.The median values of the four water balance components in the average year correspond closely to the overall yearly means presented in Table 4, confirming that the average year represents well the overall average, both in terms of rainfall characteristics and for the water balance components.The figure also shows the deterministically obtained values illustrating, again, differences between those values and the stochastic medians or means.

CONCLUSIONS
The deterministic and stochastic estimation of plant available water and the simulation of soil water balance components for a 30-year period in seven Brazilian soils cropped with maize allowed us to conclude that (1) hydraulic parameter uncertainty strongly propagates into the estimates of plant available water despite a relatively small uncertainty in the threshold water contents; (2) variations in rainfall characteristics are more determinant for bottom flux and evaporation than variations in soil properties; on the other hand, variations in soil properties are more determinant for transpiration and runoff; (3) average values of stochastic simulations for soil water balance components may differ significantly from the results of deterministic simulations; and (4) the complex interaction between soil water retention and conductivity parameters, soil water dynamics and related processes like root water uptake, evaporation, infiltration, and drainage, makes the statement of simple relations between those parameters and soil water balance components unfeasible.

A C K N O W L E D G M E N T S
The research leading to this publication was supported by the São Paulo Research Foundation (FAPESP), projects 2022/00853-2 and 2020/07294-3.Further acknowledgments to the Kirkham Foundation for covering the APC charges.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

F
Monthly distribution of rainfall and minimum and maximum temperatures observed in Piracicaba, Brazil, between 1992 and 2022.Error bars indicate ± the standard deviation.

F
Average profile water content at the limiting point (θ LP ) versus field capacity (θ FC ) estimated for the 10 4 stochastic realizations for soil A. Dashed lines are isoRAWs and numbers in italic represent the corresponding readily available water (RAW; in cm) considering a 0-60 cm profile.

F
Cumulative frequency of readily available water (RAW) and total available water (TAW) in two of the soils (soil D and G) considering a 0-60 cm profile.T A B L E 5 Simulated values of bottom flux in absolute values (|BF|, cm) and transpiration (T, cm) in the wet, average, and dry year for each of the seven soils obtained stochastically (averages of 10 4 stochastic realizations ± respective standard deviations (stoch ± σ r ) and deterministically (det.).
Location, particle size distribution, texture class, and soil order for both layers of the seven analyzed soils.
T A B L E 1 Average simulated values of the soil water balance components bottom flux (BF, cm; negative values designating a downward flow), evaporation (E, cm), transpiration (T, cm), and runoff (R, cm) considering 30 years, seven soils, and 10 4 stochastic realizations.Yearly means (mean y ) and standard deviations (σ y ) from 30 yearly mean values of 10 4 stochastic realizations.Soil means (mean s ) and standard deviations (σ s ) from the seven evaluated soils.
T A B L E 4 Transformed parameters with respective standard deviations and correlation matrix, and untransformed (VGM mean ) parameters for both layers of the seven analyzed soils.θ 2 and θ 1 correspond to pressure heads of −10 and −1000 cm, respectively.VGM, Van Genuchten-Mualem.