the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 12 Jun 2025
Soil physico-chemical indicators for ecosystem services: a focus on water regulation
Abstract. This study investigates the intricate relationship between soil properties and water-related processes, with a focus on their collective impact on ecosystem service provision, particularly water regulation. Conducted in three diverse regions Marchfeld (Austria), Bologna (North Italy) and Rmel (Tunisia), the research aims to identify key soil properties that influence water infiltration (INF), groundwater recharge (GWR), and crop water stress indexes (CWSI). Key soil characteristics such as saturated hydraulic conductivity (KS ), available water content (AWC), bulk density (BD), saturated water content (θs ), organic matter (OM), clay content and soil depth were analyzed for their role in regulating water movement and the overall hydrological balance. Pairwise correlation and multiple linear regression analyses were used to assess the interactions among soil water balance processes and soil properties. The results reveal significant variations between regions in terms of the factors that control infiltration, groundwater recharge, and CWSI. For example, in Marchfeld infiltration showed a strong positive correlation with BD (r = 0.74, p < 0.001), while CWSI had the most significant negative correlation with soil depth (r = -0.35, p < 0.001). Futhermore, multiple linear regression models were developed to assess the relevance of the different soil properties and of their interactions on the components of the soil water balance. As an example, in Marchfeld, the model for infiltration (r = 0.79, p < 0.001) was highly predictive, incorporating Clay, OM and soil depth. These results emphasize the critical role of key soil properties KS , AWC, BD, OM, clay content, θs and soil depth in controlling soil water processes. The study highlights the value of using these properties in predictive models to inform water management practices to optimize crop performance and soil conservation in different agricultural settings.
- Preprint
(1424 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 27 Oct 2025)
-
RC1: 'Comment on egusphere-2025-1927', Anonymous Referee #1, 28 Jul 2025
reply
-
AC1: 'Reply on RC1', Marialaura Bancheri, 29 Aug 2025
reply
Referee #1, 28 Jul 2025
General comments
The authors studied the link between soil physical and chemical parameters and indicators of hydrological processes and ecosystem services. They performed an extensive simulation study with the hydrological FLOWS model for 3 contrasting regions regarding soil variability and climate. They showed that the most relevant soil properties for soil functions and ecosystem services depend on the environmental conditions in a region, which is an important innovation.
The paper reads fluently, the methodological setup is sound, and the results are well presented in Tables. However, I still have questions and remarks, as described below.
A. We would like to sincerely thank the Reviewer for their careful reading of our manuscript and for their constructive comments and suggestions. Their feedback has been very helpful in improving the clarity, quality, and overall rigor of the work.
Particular comments
Q. Figure 2: For each site, a different soil legend is used. This makes comparison hard. Can you create a figure with a single uniform soil legend?
A. Although we agree with the suggestion of the reviewer, the adoption of a single legend for the three sites would result in a loss of detail in the spatial distribution of the different soil types presented in the three areas. This is due to the fact the three soil maps have a different reference scale: 1:25,000 for Marchfeld (detailed soil map, where different anthropic soils units are identified and mapped), 1:50,000 for Bologna (semi-detailed soil map, where mapping units are pedolandscapes with distinct soil series) and 1:100,000 for Rmel (reconnaissance soil map, where mapping units are in most cases associations or complexes of contrasting soil groups). This results in different levels of taxonomic details for the mapping units of the three sites, which in practice hamper the adoption of a single legend without reducing the detail of the information and the soil variability depicted in the maps.
Q. Figure 3: In the Figure of FLOWS, you show the groundwater. Does FLOWS consider capillary rise?
A. Yes, FLOWS consider capillary rise of both water and dissolved solutes. FLOWS is based on the numerical solution of the Richards equation and water fluxes are calculated according to the hydraulic gradient. In the case of a shallow water table, thus close the bottom of the simulated soil domain and interacting with the soil profile, upward water fluxes (capillary rise) will develop anytime the hydraulic head at the water table will be higher than that at the simulation node just above the water table.
Q. Also related to this Figure: “lateral drainage” is drawn in the unsaturated zone. As the majority of lateral drainage occurs in the saturated zone, better to depict “lateral drainage” there.
A. There appears to have been a misunderstanding regarding the use of the term drainage. In the FLOWS model, “drainage” does not refer to lateral drainage processes naturally occurring anytime there is a lateral hydraulic gradient. Actually, FLOWS is a 1D model and would not allow for calculating these lateral fluxes. Instead, in FLOWS this term has been used to identify fluxes to artificial drains, calculated according to the Hooghoudt’s theory, which are included in the Richards equation as a sink term.
For clarity, in the figure the term drainage was changed to artificial drainage, to distinguish it from the vertical drainage and to avoid any confusion with lateral Darcian flow.
Q. Line 130: The FLOWS model plays an important role in your analysis. Can you mention studies that showed the reliability of FLOWS?
A. The model has been used for several years in several papers related to water and solute balance, even in complex soil profiles, both at the theoretical (see for example, Coppola et al., 2009; 2012; 2015a) and applicative levels (see for example Coppola et al., 2015b; 2019a, 2019b; Hassan et al., 2024; Wang et al., 2013; Porru et al., 2024).
Coppola, A., M. Abdallah, G. Dragonetti, P. Zdruli, and N. Lamaddalena. 2019a. “Monitoring and Modelling the Hydrological Behaviour of a Reclaimed Wadi Basin in Egypt.” Ecohydrology 12, no. 4: e2084. https://doi. org/ 10. 1002/ eco. 2084.
Coppola, A., G. Dragonetti, A. Sengouga, et al. 2019b. “Identifying Optimal Irrigation Water Needs at District Scale by Using a Physically Based Agro- Hydrological Model.” Water 11, no. 4: 841. https:// doi. org/10. 3390/ w1104 0841.
Hassan, S.B.M., Dragonetti, G., Comegna, A., Lamaddalena N., Coppola A, 2024. Analyzing the role of soil and vegetation spatial variability in modelling hydrological processes for irrigation optimization at large scale. Irrig Sci (2024). https://doi.org/10.1007/s00271-023-00882-7
Coppola, A., Basile A., Comegna A., Lamaddalena N., 2009. Monte Carlo analysis of field water flow comparing uni- and bimodal effective hydraulic parameters for structured soil, Journal of Contaminant Hydrology (2008), doi:10.1016/j.jconhyd.2008.09.007.
Coppola, H. H. Gerke, A. Comegna, A. Basile, V. Comegna, 2012. Dual-permeability model for flow in shrinking soil with dominant horizontal deformation. Water Resources Research, Vol. 48, W08527, doi:10.1029/2011WR011376.
Coppola A., A. Comegna, G. Dragonetti, H. H. Gerke, A. Basile, 2015a. Simulated Preferential Water Flow and Solute Transport in Shrinking Soils. Vadose Zone J. doi:10.2136/vzj2015.02.0021
Wang X., G. Quan, Y. Pan, R. H.,Y. Zhang, A. Tedeschi, A. Basile, A. Comegna, A. Coppola, R. de Mascellis, 2013. Comparison of hydraulic behaviour of unvegetated and vegetation-stabilized sand dunes in arid desert ecosystems. ECOHYDROLOGY (ISSN:1936-0584). DOI: 10.1002/eco.1265;
Coppola A., N. Chaali, G. Dragonetti, N. Lamaddalena and A. Comegna, 2015b. Root uptake under non-uniform root-zone salinity. Ecohydrol. (2014). Wiley. DOI: 10.1002/eco.1594
Porru M.C., Hassan S. B. M., Abdelmaqsoud M., Vacca A., Da Pelo S., Coppola A., 2024. Using index and physically-based models to evaluate the intrinsic groundwater vulnerability to non-point source pollutants in an agricultural area in Sardinia (Italy). Front. Water. 6:1399170. doi: 10.3389/frwa.2024.1399170
In addition, two manuscripts specifically addressing an inter-code comparison with the HYDRUS model are in preparation and will be submitted shortly to Computational Geosciences. Their provisional titles are:
FLOWS: A model to simulate water and solutes fluxes in agricultural and environmental systems. 1. Theory and applications with inter-code benchmarking
FLOWS: A model to simulate water and solutes fluxes in agricultural and environmental systems. 2. Modules for specific applications in agro-hydrological modeling
Q. Line 147: At which soil depth did you apply the free drainage option?
A. The free drainage boundary condition was applied at the bottom of each soil profile, with the depth corresponding to the specific profile considered in the simulations. We will clarify this point in the revised manuscript.
Q. Line 159: You define here the potential soil evaporation, which in general is much larger than the actual soil evaporation, an important part of your water balance. How did you calculate the actual soil evaporation?
Q. Line 165: Is potential root uptake also assumed to be proportional to root density?
A. The two questions are somehow related. FLOWS uses the reference evapotranspiration and a crop coefficient to calculate the potential evapotranspiration, ETp. This, in turn, is separated in the potential evaporation, Ep, and potential transpiration, Tp, components by using the Beer’s law, which includes an attenuation factor for radiation penetration through the vegetation.
As for actual evaporation, Ea, FLOWS calculates the maximum upward flux at the soil surface, Emax, that the soil conditions may allow, based on the hydraulic gradient between the top boundary and the first node in the soil profile, as well as on the hydraulic conductivity in that numerical compartment. Accordingly, the code uses the following condition for calculating Ea: if Ep<=Emax, Ea=Ep; if Ep>Emax, Ea=Emax.
As for actual transpiration, Ta, FLOWS uses a macroscopic root uptake (see for example, Feddes, Kowalik, and Zaradny 1978; Feddes and Raats 2004). This approach distributes the potential transpiration, Tp, along the different simulation (numerical) nodes within the rooted depth (which is an input for the model) according to a root density distribution function, g(z). The user may opt for different root distributions, such as uniform, triangular, logistic, etc. The fraction of Tp in a node multiplied by the node compartment depth, delta_z, provides the potential sink, Sp, at that node. In other words, the integral of Sp over the rooted depth corresponds to Tp.
To obtain Ta, FLOWS firstly calculates an actual sink term, Sa, by multiplying the Sp by an uptake reduction function, α(h,hos). The sink reduction factor accounts for the possible water stress (related to soil water pressure head, h) and osmotic stress (related to the osmotic potential, hos) experienced by roots at different simulation nodes within the root zone. As the reduction uptake factor is related to the local (the specific numerical node) conditions within the root zone, there may be nodes where roots are drawing water at the potential conditions (α=1 and Sa=Sp) and stressed nodes (α<1 and Sa<Sp). The integral of Sa over the rooted depth corresponds to Ta, which may thus be <= to Tp, depending on the stress conditions along the root zone.
Accordingly, details on how the model calculates actual evaporation and transpiration will be provided in the revised text.
Q. Lines 174-183: You select for each site a median year based on the amount of rainfall. Did you also check whether the rainfall distribution over the year was representative of the site's climate?
A. In addition to selecting the median rainfall year for each site, we also examined the intra-annual rainfall distribution to ensure that the selected year was representative of the typical seasonal rainfall pattern of each site.
Q. Precipitation input: At which interval did you provide input of rainfall fluxes? In order to simulate surface runoff, daily rainfall intensities are generally not suitable.
A. The model simulations were carried out at a daily time step, with precipitation provided as daily amounts. We acknowledge that this temporal resolution does not allow for the explicit simulation of surface runoff processes driven by short-duration, high-intensity rainfall events. However, the purpose of our work was to simulate the soil–water balance at the daily scale, focusing on infiltration, evapotranspiration, and soil moisture dynamics, rather than a single event-based runoff generation. Moreover, precipitation data at shorter time steps are rarely available at the large spatial scale of our three case studies, which further justifies the use of daily inputs.
Q. Lines 184 and 185: The data of soil texture and organic matter: to which soil layer do they belong? Did you assume these data to be uniform over the entire root zone depth (max 80 cm)?
A. The data on texture, bulk density, and organic matter content were available for each horizon of the investigated soil profiles. For each horizon, the hydraulic parameters were then derived by applying the PTFs. Finally, for the purpose of the statistical analysis, we computed (i) the weighted mean of these properties over the horizon depths, and (ii) the arithmetic mean across all soil profiles within each site, in order to obtain overall spatial statistics.
Q. Line 187: The crop coefficient is defined depending on the stage of crop development. How did you derive these stages?
A. The crop development stages were determined by integrating the FAO56 guidelines with expert knowledge, literature, and site-specific field observations. This approach ensured a robust and site-appropriate parameterization of crop development for each case study. For example, based on the typical crop calendar and the growth stages described in FAO56, flowering occurs earlier at RMEL than at Marchfeld and Bologna; while sowing takes place earlier at Marchfeld and Bologna compared to the Rmel site. To further strengthen transparency, we will provide the representative development stages for each site in the supplementary material and revise the main text accordingly.
Q. Table 2: The water balances of Bologna and Rmel do more or less close. However, the fields in Marchfeld, on average gained 94 mm of water. This is a substantial amount compared to other water balance components and a bit strange for a median year. Please clarify this in the text.
A. From Table 2, the higher average water storage (91.8 mm) observed in Marchfeld soils compared to RMEL (14.4 mm) and Bologna (-3.3 mm) can be attributed to their greater soil depth (Table 1) and the resulting right-skewed distribution of soil profile depths, which allows more infiltrated water to remain within the profile. Moreover, the use of a single climatic year prevents proper closure of the water balance as would be expected when extending the simulations over longer periods.
Q. Line 274: Minimum and maximum values are not depicted in Table 2.
A. We have now updated the table to include the minimum and maximum values for each variable, providing a more complete picture of the variability across sites.
Q. Lines 417 – 433: Important that you mention the ambiguity on the role of the saturated conductivity. You might add here that the simulated amount of runoff not only depends on the saturated conductivity but also on the time intervals with which rainfall is provided. In general, daily intervals underestimate the runoff amounts.
A. Thank you for your comment. We will add the following text to clarify this aspect:
“This analysis was conducted at a daily scale; however, for certain terms of the water balance, such as potential runoff production, a finer temporal resolution would be more appropriate. This is particularly relevant in the context of increasing high-intensity, short-duration rainfall events associated with climate change (Prein et al., 2017).”
Prein, A. F., Rasmussen, R. M., Ikeda, K., Liu, C., Clark, M. P., & Holland, G. J. (2017). The future intensification of hourly precipitation extremes. Nature climate change, 7(1), 48-52.
4 Discussion: Mention important limitations of your study:
- Q. Your hydrological results are based on median weather years. Commonly, in these years, hydrology functions well. The problems might arise in extreme dry or wet weather years. A thorough analysis of soil functions and ecosystem services should include more extreme weather years.
A. We fully agree that assessing soil functions and ecosystem services under extreme weather conditions is essential. While the current analysis is based on median weather years to provide a baseline characterization, we are already working on extending the assessment to include both dry and wet extremes. This will allow us to better capture the variability and potential vulnerabilities of soil functions under climate variability and change.
Accordingly, the following text will be added to the manuscript in the discussion section at line 455:
“It should be noted that, in assessing the complex interplay between meteorological variables and soil properties, this study was conducted using a median year, which provides baseline information on climate–soil interactions. Clearly, extending such analyses over longer time spans, under more extreme climatic conditions, and at finer (e.g., hourly) temporal resolution would provide additional and complementary insights.”
- Q: You focus on the durum wheat crop. In the case of other crops or multiple crops per year, the results will be different.
A: Certainly, the results would vary with different crops or multiple cropping systems. We focused on durum wheat because it is a widely cultivated rainfed crop in the three study areas, and its growth cycle is well-documented. This allowed us to rely on robust phenological information and expert knowledge to accurately define its development stages. Moreover, the absence of irrigation simplifies the analysis by reducing variability due to water management, making it a suitable reference case for our assessment framework.
- Q: Soil layering, which is explicitly mentioned in the introduction, is not included in the analysis.
A: Soil layering was in fact considered in the analysis, but we acknowledge that this was not sufficiently explained in the manuscript. As also noted in the previous comment, we will revise the text to clearly indicate how soil layering was incorporated into the model setup and analysis.
Detailed comments
Q: Lines 110 – 114: Current text is wrongly phrased and repetitive. Define the climate of site Rmel more concisely and clearly.
A: Corrected
Q: Line 134: “Genucthen” should be “Genuchten”.
A: Corrected
Q: Line 138: The symbol “h” was already defined in line 134. Define each variable once, upon first usage in the text.
A: You’re generally right. We are going to revise and delete the repetition of the variables when it is possible, according to the specific contexts.
Q: Lines 189 and 190: Include references.
A: Corrected
Q: Line 191: Include units.
A: Corrected
Q: Line 134: “groupf” should be “group”.
A: Corrected
Q: Table 2: In the caption, clarify CWSI.
A: Corrected
Q: Table 3: In the caption, clarify SD.
A: Corrected
Citation: https://doi.org/10.5194/egusphere-2025-1927-AC1
-
AC1: 'Reply on RC1', Marialaura Bancheri, 29 Aug 2025
reply
-
RC2: 'Comment on egusphere-2025-1927', Anonymous Referee #2, 27 Aug 2025
reply
General remarks:
The authors present a study on the use of the FLOWS model for identifying the correlation between soil indicators and ecosystem services. The manuscript is properly written, describing how the study is performed adequately. However, there are a number of issues that before this study can be published.
- At lines 39 – 43 you state that the current literature is in a simplified manner showing integration of soil indicators and their overall effect on soil functions, soil processes and ecosystem services. This is a rather simplified summary of the literature review, since there are many studies actually showing an integrative, holistic approach of soil indicators (e.g. Drobnik et al., 2018, Schoenholtz et al., 2000, Dominati et al. ,2010, Dominati et al., 2016, Constantini et al., 2016, to mention a few). In fact, you are not extending your analysis to a holistic approach. You use a limited set of physical soil indicators (soil hydrological parameters, partly even taken from PTFs form soil texture data) to explain their effects on a limited set of ecosystem services (groundwater recharge, infiltration, and the calculated rest products of these). What is your contribution to your statement that current studies show a simplified manner in integrating soil indicators and ecosystem services? Clarify this in your introduction and results/discussion accordingly.
- On the methodology part: you have used the FLOWS model to calculate your identified ecosystem services (RUN, INF, GWR, CWSI). Basically, this study is a sensitivity analysis of the model, without a proper description of the model itself. You indicate that you use a ‘representative median year’ for running the model (line 170), but not what the calculation timesteps of your runs are. Coppola et al.(2019) used rain events for calibrating/validating the model, with the, among others, Ks as calibration factor. You state at line 353 that Ks has a moderate to weak correlation for infiltration (and even more pronounced at line 418). But in what way is this due to the configuration of your model-design? How does the timestep affect this, or the number of profile layers in the model setup? Is there a top-layer and a more compact bottom layer in your profile? In addition, the Ks you use is not the actual Ks, since macropore flow is not taken account for by the MvG equations. There is a limited number of publications on the use of the FLOWS model, so more information is essential to understand its concept. Moreover, you do not calibrate or validate your model runs, so it remains unclear if the magnitude of the results is realistic. Explain your model configuration, validation and perform a sensitivity analysis.
Other remarks:
Line 68: you are not evaluating chemical soil parameters, only physical. OM is often considered both chemical and physical, but still leaves you with only one (partly) chemical.
Line 75: refer to the original paper of the model (Coppola et al., 2019), although in that paper the model is also not properly described. This makes understanding the model difficult.
Line 111 – 112: Both rainfall and precipitation are mentioned, use the same phrasing.
Line 134: Typo in Genuchten
Line 199 (and elsewhere): you state at places that the runoff is calculated. But the model does not calculate runoff. Actually, at line 200 you mention the correct terminology: ‘potential runoff’. Basically, this is the outcome of the equation:
Potential runoff = Precipitation – Interception – infiltration.
The remainder is ponding, potentially available for runoff (but also for evaporation or infiltration in the next timestep of calculations). Change this term accordingly.
Line 207: Typo in ‘group’
In Results: Section 3.1 should go to ‘methods’, it describes the study area and input data.
Table 1: You present the MvG parameters, as a result of the PTFs given by Wôsten et al., 2001. I assume you are aware of the very low R2 of most of those values, making your study even more difficult without having the option for validating your model results. Describe this in your discussion, how does this fit affect your results.
Line 289: The discrepancy is not due to runoff, since runoff (better: ponding) is the final component of the calculations. So, your statement that runoff removes substantial portion of the precipitation is not true, since infiltration is calculated first.
Citation: https://doi.org/10.5194/egusphere-2025-1927-RC2 -
AC2: 'Reply on RC2', Marialaura Bancheri, 10 Sep 2025
reply
Referee #2, 27 Aug 2025
General remarks:
The authors present a study on the use of the FLOWS model for identifying the correlation between soil indicators and ecosystem services. The manuscript is properly written, describing how the study is performed adequately. However, there are a number of issues that before this study can be published.
A. The authors sincerely thank the thoughtful and constructive feedback on our manuscript. They appreciate the recognition of the overall quality of the work and the adequacy of our study description. The authors will carefully address the issues raised by the Reviewer to further improve the clarity, rigor, and contribution of our study before resubmission.
Q. At lines 39 – 43 you state that the current literature is in a simplified manner showing integration of soil indicators and their overall effect on soil functions, soil processes and ecosystem services. This is a rather simplified summary of the literature review, since there are many studies actually showing an integrative, holistic approach of soil indicators (e.g. Drobnik et al., 2018, Schoenholtz et al., 2000, Dominati et al. ,2010, Dominati et al., 2016, Constantini et al., 2016, to mention a few). In fact, you are not extending your analysis to a holistic approach. You use a limited set of physical soil indicators (soil hydrological parameters, partly even taken from PTFs form soil texture data) to explain their effects on a limited set of ecosystem services (groundwater recharge, infiltration, and the calculated rest products of these). What is your contribution to your statement that current studies show a simplified manner in integrating soil indicators and ecosystem services? Clarify this in your introduction and results/discussion accordingly.
A. At lines 39-43 we wanted to remark that the inherent complexity of soil functions resulting in the processes underpinning the delivery of soil-based ecosystem services hamper to some extent the efficacy of simplified indicators which, by definition, are static being based on a limited set of soil properties. If from the one hand holistic approaches of soil indicators have allowed the definition of theoretical frameworks for the integration of soil-based ecosystem services into policy and land planning, and actually being quite successful at doing so, from the other there is the risk to overlook the actual complexity of field dynamics and their temporal variability. Our contribution was not intended to extend the analysis of modelling results in the three case study areas to a more holistic approach, but rather to provide insight into the relationships between the outcomes of a dynamic process-based model and a set of soil properties frequently used as proxies for the assessment of regulation ecosystem services, focusing on water regulation in different climatic conditions.
We have clarified the aim of our contribution by modifying the text in the Introduction and Results/Discussion sections, as suggested by the Reviewer.
Q. On the methodology part: you have used the FLOWS model to calculate your identified ecosystem services (RUN, INF, GWR, CWSI). Basically, this study is a sensitivity analysis of the model, without a proper description of the model itself. You indicate that you use a ‘representative median year’ for running the model (line 170), but not what the calculation timesteps of your runs are. Coppola et al. (2019) used rain events for calibrating/validating the model, with the, among others, Ks as calibration factor. You state at line 353 that Ks has a moderate to weak correlation for infiltration (and even more pronounced at line 418). But in what way is this due to the configuration of your model-design? How does the timestep affect this, or the number of profile layers in the model setup? Is there a top-layer and a more compact bottom layer in your profile? In addition, the Ks you use is not the actual Ks, since macropore flow is not taken account for by the MvG equations. There is a limited number of publications on the use of the FLOWS model, so more information is essential to understand its concept. Moreover, you do not calibrate or validate your model runs, so it remains unclear if the magnitude of the results is realistic. Explain your model configuration, validation and perform a sensitivity analysis.
A: This comment from Reviewer 2 covers several topics. For clarity, we have structured our reply into two parts: the first addresses the overarching questions, while the second focuses on specific points.
A. Overarching questions:
On the comment regarding the ‘representative median year’ used for running the model (line 170), as already replied to Reviewer 1, in addition to selecting the median rainfall year for each site, we also examined the intra-annual rainfall distribution to ensure that the selected year was representative of the typical seasonal rainfall pattern of each site. As described in Section 2.3, the representative median climatic year was calculated for each of the three sites. The selected representative median year corresponds to the year within the full observation period whose annual precipitation was closest to the median of the long-term dataset at each site. This approach ensures that the chosen year reflects typical precipitation and evapotranspiration conditions, making it representative of the site’s median hydro-climatic regime.
One of the reviewer’s requests concerns additional details on the model. While it is not feasible to include all equations and processes implemented in the model within the manuscript, we will clarify several aspects that were identified as unclear in Section 2.2 and provide a link to the model manual for further technical details (Coppola et al., 2025). In addition, as already noted in our reply to Reviewer 1, two manuscripts including inter-code comparisons with HYDRUS-1D are currently in preparation and will soon be submitted.
Coppola, A., Hassan, S. B. M., & Pili, F. (2025). FLOWS - FLOw of Water and Solutes in agricultural and environmental systems (WEB Version 1.0). Zenodo. https://doi.org/10.5281/zenodo.17076523
The second general comment concerns what we consider a fundamental rationale for this work, which may not have been adequately and convincingly conveyed in the initial manuscript. This study is not conceived as an exercise in applied modelling in the strict sense. Rather, its primary purpose is to emphasize—particularly for the community of soil scientists focused on soil health and ecosystem services, but less familiar with soil hydrology—that relationships between indicators and hydrological processes are inherently complex and indirect. Such relationships are influenced by multiple interacting factors (e.g., soil layering, upper and lower boundary conditions) that cannot be captured by a simple one-to-one indicator–process assessment.
Finally, we do not consider this work reducible to a mere sensitivity analysis of model parameters, mainly because our evaluation is not based on perturbations of direct model inputs (such as van Genuchten parameters or LAI). Instead, it relies on widely used indicators of water regulation ES, which may coincide with direct model parameters (e.g., Ksat) or may not (e.g., AWC).
- Specific comments: (on the sentence Ks shows only moderate to weak correlation with infiltration)
- Time steps used in the simulation runs: Information on the daily time step used in the simulations is already reported in Section 2.2.2 Boundary and initial conditions for water flow. To improve clarity, we will restate this point when presenting the simulation results. In addition, as also discussed in our response to Reviewer 1, we acknowledge that for specific terms such as runoff generation, a finer (e.g., hourly) time step would be more appropriate, particularly in the context of high-intensity, short-duration rainfall events.
- On the role of soil layering: Each soil profile was discretized into multiple horizons based on available site data. This stratification is a key factor influencing infiltration, as water redistribution across horizons reduces the direct influence of topsoil Ks on overall infiltration. The soil layering represents a clear example of the soil complexity highlighted in the manuscript and the need of a modelling approach. For instance, a sand horizon overlying clay produces completely different infiltration and redistribution dynamics compared to the reverse arrangement, with outcomes that cannot be predicted a priori. For this reason, a physically-based modelling approach can adequately capture such interactions and provide reasonable estimates of water balance components. Conversely, identical soil properties arranged in a different layering sequence, or under a different climate, may lead to markedly different results.
Regarding the Reviewer’s question on the presence of a top layer overlying a more compact bottom layer, given the large number of soil profiles considered, we observed a variety of situations. Not all soils follow the classical sequence of increasing compaction with depth. Instead, below the ploughed Ap horizon, several horizons with different degrees of compaction are present, resulting in a range of stratifications across the profiles.
Finally, since this information on the number of soil horizons was missing from the original text, we have now added the following sentence in ‘Section 2.3 - Input dataset for model simulations’: “The number of horizons reported in the soil maps of the three areas ranges from a minimum of 1 to a maximum of 5. On average, Marchfeld, Bologna and Rmel profiles include 3 horizons.
- Ks influenced by model design/configuration: Certainly, the correlation between Ks and the different processes (i.e., infiltration) is influenced by the model design/configuration. For instance, a bucket-type model with a homogeneous soil profile would yield markedly different results compared to a physically based approach, as used in our study, where processes are simulated with state-of-the-art numerical methods. In particular, the physically-based approach with the Richards' equation explicitly accounts for real driving force, soil layering, boundary conditions, flux-to-runoff, and the dynamic redistribution of water within the profile, thereby reducing the apparent direct control of topsoil Ks alone and providing more realistic insights into infiltration and water balance dynamics.
- Ks influenced by macropores: The reviewer is certainly correct that the Ks used in our simulations does not represent the field-saturated conductivity (Ksf ). This limitation is well recognized by the authors, and for this reason we deliberately denoted it simply as Ks , in line with common practice to distinguish between the two. Furthermore, considering that our analysis involved simulations across 315 soil profiles, it is practically impossible—at least to our knowledge—to obtain consistent information on macropore presence and their parameterization at such a scale. Therefore, our focus was placed on matrix conductivity, while acknowledging that macropore flow is an important process at the local scale but remains beyond the scope and feasibility of this study.
- Implications for the realism of the magnitude of results: Explain your model configuration, validation and perform a sensitivity analysis. We appreciate the Reviewer’s concern regarding the validity of the model and the request for further information. To address this, we note that the full manual—containing all equations and details of the model design—together with the scientific papers in which the model has already been published and validated, provides a solid basis for transparency and reproducibility. In addition, Chapter 5 of the manual (Coppola et al., 2025) evaluates the reliability of FLOWS by comparing its results with those obtained from HYDRUS-1D simulations, further supporting confidence in the model’s performance.
With respect to the request for a sensitivity analysis of the model parameters, we fully agree that this is a valuable exercise. Such an analysis is indeed in preparation and will further strengthen the model’s applicability.
Other remarks:
Q: Line 68: you are not evaluating chemical soil parameters, only physical. OM is often considered both chemical and physical, but still leaves you with only one (partly) chemical.
A: We agree that our study mainly focuses on soil physical parameters (e.g., texture, bulk density), while organic matter (OM) can be considered both a chemical and a physical property. To avoid ambiguity, we revised the text accordingly and now refer to “soil parameters” instead of “chemical and physical parameters.”
Q: Line 75: refer to the original paper of the model (Coppola et al., 2019), although in that paper the model is also not properly described. This makes understanding the model difficult.
A: See our answer in the “Overarching questions” section
Q: Line 111 – 112: Both rainfall and precipitation are mentioned, use the same phrasing.
A: We will use the term rainfall throughout this work.
Q: Line 134: Typo in Genuchten
A: Thanks, we will correct the Genuchten in the text
Q: Line 199 (and elsewhere): you state at places that the runoff is calculated. But the model does not calculate runoff. Actually, at line 200 you mention the correct terminology: ‘potential runoff’. Basically, this is the outcome of the equation:
Potential runoff = Precipitation – Interception – infiltration.
The remainder is ponding, potentially available for runoff (but also for evaporation or infiltration in the next timestep of calculations). Change this term accordingly.
A. The model does not compute runoff as a residual of the water balance; rather, it explicitly simulates both Hortonian and Dunnian runoff generation processes. We will clarify this point in the model description. To prevent misunderstanding, we have replaced the term “runoff” with “flux-to-runoff,” which more accurately represents the model’s approach. FLOWS is a one-dimensional model and does not formally simulate routing runoff. Furthermore, the model allows the user to specify a maximum acceptable ponding depth (e.g., 10 cm): any flux exceeding this threshold is considered flux-to-runoff. After the rainfall ceases, the ponded water gradually infiltrates into the soil profile.
Q: Line 207: Typo in ‘group’
A: Thanks, we will correct the group in the text
Q: In Results: Section 3.1 should go to ‘methods’, it describes the study area and input data.
A: In Subsection 3.1, as indicated in the introduction to the results, we aimed to present the soil data used in the simulations, emphasizing the study area’s diverse pedological characteristics. For this reason, we considered it appropriate to include this subsection in the results. At the same time, we recognize that the table itself is more closely related to the description of the data and procedures, and therefore we have moved the table to Section 2 (Materials and Methods), while keeping the interpretative text in the results. This way, both clarity and consistency with scientific conventions are maintained.
Q: Table 1: You present the MvG parameters, as a result of the PTFs given by Wôsten et al., 2001. I assume you are aware of the very low R2 of most of those values, making your study even more difficult without having the option for validating your model results. Describe this in your discussion, how does this fit affect your results.
A: The reviewer is certainly correct in pointing out the limitations of relying on PTFs. This remains a fundamental challenge, especially when addressing problems at a regional rather than a local scale. As we have already written in our response to Reviewer 1, we will strengthen the discussion section to more explicitly acknowledge and elaborate on the limitations associated with the use of PTFs, so that the reader is fully aware of their implications in the present context.
Q: Line 289: The discrepancy is not due to runoff, since runoff (better: ponding) is the final component of the calculations. So, your statement that runoff removes substantial portion of the precipitation is not true, since infiltration is calculated first.
A: This question appears to stem from a misunderstanding of how the model calculates runoff (now referred to as flux-to-runoff). We kindly refer the Reviewer to our previous response on this issue, where the procedure is explained in detail.
Citation: https://doi.org/10.5194/egusphere-2025-1927-AC2
-
RC3: 'Comment on egusphere-2025-1927', Anonymous Referee #3, 20 Sep 2025
reply
The manuscript presents and elaborate and well structured study in a nicely written and structured way, for which I congratulate the authors. However, I do think that there several key issues that would need to be carefully considered before this submission can be published. Some of the comments call for potentially substantial amount of additional work and re-analysis of the data.
- On the description of the modeling case: I somewhat lack a proper description of the modeling process, parameterization, etc.: a clear delineation and discussion of which of the reasonably important processes are (or are not) modelled by the model. Was the full functionality of the model used or did data limitation or some other reason reduce what processes were actually considered? For instance, preferential flow is a key process to consider given the aim of the paper. I also generally lack clarity on the parameterization of the model, how many profiles have been modeled, as well as any mention of validation of the model for the purpose or location(s).
- On the matter of scale: To me it feels that there is a lack of connection established between the scale of modeling and the scale at which conclusions are made: The model is apparently a 1D simulation model that has been run on individual soil profiles. However, “spatial differences” have been referred to (around L255), while there has not been any information provided on how the upscaling to 100s of km2 has been handled. Each area is obviously heterogeneous in terms of soil type and depth. I think this should be discussed, along with the limitations of the chosen approach.
- On using only the average year in the simulations: I find working with just a single “average” year to be a weakness of the study. Most of the events of concern – be that runoff, erosion, pollution, drought, etc. – are usually the products of either extreme events or extreme seasons/years. This is unfortunately not considered in the study, and I feel that the given advice is not based on the best available information. Working with a distribution of weather conditions (years) at each site would also allow the deduction of correlations between various factors in a much deeper way, allowing the study’s conclusions to be much better supported. Given that the data are available, the workflow would just need to be extended, and the results analysed at a higher dimensionality. I think the best value for the reader would be if the analysis was rerun on all available years and the correlation analysis was deeper, accounting for interactions and distributions that look into e.g. precipitation distribution (same soil under wetter and drier conditions) and soil depth as well, which is crucial when capacity-driven metrics are of concern.
- On the diversity of soil depths and their handling: Soil depth is obviously one of the keys to the findings, whereas it has not really been addressed as such. It has mentions, but it is a key factor that I believe requires much more attention. The study mentions e.g. drainage and ground water recharge. At what depth are they defined in this study? How are these factors handled when the soil is shallow?
- On some of the obvious uncertainties that are not being acknowledged and handled: Let me pick an example. I think the conclusions made regarding Ksat vs. infiltration at L417-429 are not thought through sufficiently and therefore are a bit careless. Essentially the issue is a set of limitations that originate from possibly several sources. In this study Ksat has been estimated by a generic PTF, which smoothens differences and is generally error prone. The other variable, infiltration is a result of modeling, which is prone to its own set of uncertainties and errors. Neither of them are measured here. Ksat and infiltration are not the culprits, but the lack of our ability to account for them is. Such sets of uncertainties should be acknowledged and discussed, and the findings placed in context, and the conclusions handled accordingly.
- On advancing our understanding: The authors cite that past studies generally did not go to sufficient depth (L40-41), however, as a reader I do not really see how this study advances the state-of-the-art all that much. Much of the conclusions are rather common sense and reinforce relatively basic, textbook type knowledge that is relatable to soil depth and the “average” soil’s physical and hydrological characteristics. To refer to some examples, several very logical connections can be made in Table 2 already: infiltration is strongly influenced by Ksat (and likely soil thickness), runoff is strongly connectable with prevailing precipitation intensity and soil texture type, drainage by precipitation characteristics and Ksat. If there is no rain, there is no infiltration or runoff. If the soil is more course textured, there is good infiltration and no runoff. Examples are in e.g. lines 285-286, 291-318 and 336-337.
The quantitative confirmation by modeling is appreciated of course, but I feel that the reader does not receive sufficiently impactful new information that is of interest to THEM.
The authors should generally recognize and discuss the study’s limitations that may root from the model itself, the choices made in parameterization, the scale problems, simplifications (depth, climate “average”), and also data sources. I think the value of the reported study could be enhanced a lot by both addressing several of these with improved analysis, and discussing the rest.
Some detailed comments:
L40-41: This statement invites that the authors make a very clear statement later on how their study improves on this.
L43-48: Here I recommend consulting e.g. https://doi.org/10.1016/j.earscirev.2021.103873 and https://www.nature.com/articles/s41586-018-0463-x and citing them if they are found relevant and fitting.
L151-152: Please provide more on this – how the rule/decision was made which condition to apply. Also: was there any warm-up period considered for the model?
L155-156: How were those parameters obtained?
L167: How were the h1-h4 thresholds handled? Crop specific? Fixed? Please provide info.
L170: The choice of median year largely dismisses the significance of infiltration capacity in really wet years. Please acknowledge this appropriately.
L185-186: Please use the original 1999 Geoderma reference for the HYPRES PTFs.
L207-210: At least some of the outputs and inputs/indicators are rather closely correlated. Could you please justify their handling? Example: infiltration vs. Ksat. Also: define AWC.
L231-234: Please clarify the distinction between statistical and multivariate “relationships”. Latter is also a statistical method.
L244: Is soil depth accounted for in the modeling approach? Please describe. Table 1 suggests that there have been huge differences in soil depth _within_ each study area. Soil depth of course matters a lot when it comes to the soil’s buffering capacity in e.g. infiltration capacity, total pore volume in the profile, total available water for the crop, etc.
L264-265: These differences in coarse and fine textured soil Ks values are barely different. Usually the variability of local measurements would not even make such values statistically different. I recommend double checking and also noting this in the manuscript. PTFs tend to smoothen within the available data range, and reduce such differences normally seen in the field.
L272: “all” soil profiles? Please clarify what this really means. I also note that averaging over a huge diversity of soils within each region may simply mask the hot-spots, that matter a lot in reality.
L336-337: re: stress correlations: I think these were logically expected. Bologna has the deepest soils on average, as well as the most clayey ones, and the highest amount of precipitation. Increased clay just adds to retention, and it is likely that the “h1-h2” (wet) stress limit kicks in more often, so the reduction of clay has a positive effect. In Rmel the conditions are exactly the opposite: dry region, highest sand content, shallowest soils “on average”. Then the “h3-h4” (drought) stress limit is likely often met, so any clay present will help retain water and make it available to plants. To me this feels rather simple “on average”.
L336-339: Why not run all years and then plot and correlate precipitation, clay content and CWSI? A cloud of points and deeper analysis of the potentially available data would lend much more support to the claims made.
L340: This is against classic textbook knowledge, although it has been reported by others before. The authors could attempt to go deeper in explaining it. A starting point can be https://doi.org/10.2136/sssaj2004.0055 or https://doi.org/10.2136/vzj2016.03.0021, but there are a few other references available since.
L385-388: I think the main point is not to increase some correlation coefficients, but to help understand how variable interactions help explain why certain correlations between variables and indicators are seen and how they interact, especially in “non-average” cases. I suggest trying to focus on what can be generalized for the benefit of the independent reader, who do not associate directly with the cases of Bologna or Marchfield or Rmel.
L392: What T values are referred to? Introduce the statistic used and show the data that you are referring to.
L399-401: Without depth considered, basically the same thing has been stated already on L336-337.
L400-402: On L335-336 it has been stated that in Bologna moderate correlation was found and higher clay content has been associated with increased stress. Here it is explained that only depth (i.e. not clay) showed correlation with stress. So which one is it? Running different analyses should not be the explanation. Then one of them is the incorrect method.
L416: The variability within each of the 3 areas has been greatly masked, whereas that is what would inform the reader how to judge their own cases. I don’t really see deep enough support for that in this study.
L417-429 (and also to some extent L430-439): Ksat is being singled out here, so I come to its defense to some extent. The authors should deepen their arguments at least where prior knowledge exists. Ksat of course is known to show great spatial variability that is very difficult to capture by PTFs that use only basic physical soil properties that have very limited relationship with soil structure. The estimates by PTFs are known to be smoothened (often a lot) compared to the original data support, and be often weak descriptors of the original data, no matter what is claimed about them. I ask the authors to look up the R2 value of e.g. the used HYPRES European PTFs – they will be surprised how low it is (around 0.2, as I remember). The mean Ksat values presented on L265-266 show such sign of smoothing. The problem is by far not only with Ksat and infiltration not being strong-enough correlated, but that in this study a PTF was used – and a rather generic one. Essentially, this study did not have Ksat to work with. It had some estimations of Ksat to work with – i.e. not necessarily the correct data. Second: I feel that the model’s imperfections, simplifications, simulated processes have not been thoroughly accounted for and discussed. The authors try to draw conclusions between an estimated variable and an indicator whose simulation is dependent on a multitude of imperfectly modeled processes, under “average” climatic conditions, and smeared over varied soil depths/types. I think all such uncertainties should be acknowledged, and at least some that are feasible to handle in a more detailed way should be imroved. Relevant conclusions should of course be modified (softened?) accordingly.
L440-455: I feel that the authors could have done a deeper analysis to support this section with more results, and make it more concrete. Rainfall distribution could have been characterized, the frequency of intense (defined as…?) rain events quantified, etc.
456: This is exactly what is needed – basing the study on appropriate data. I don’t want to say that a study based on PTF estimates is illegitimate, but the conclusions drawn should account for the source of errors that originate from PTF use, model use, spatial and temporal simplifications and averaging, etc.
Citation: https://doi.org/10.5194/egusphere-2025-1927-RC3 -
AC3: 'Reply on RC3', Marialaura Bancheri, 13 Oct 2025
reply
Q: The manuscript presents and elaborate and well structured study in a nicely written and structured way, for which I congratulate the authors. However, I do think that there several key issues that would need to be carefully considered before this submission can be published. Some of the comments call for potentially substantial amount of additional work and re-analysis of the data.
A: The Authors are thankful for Reviewer’s thorough and thoughtful comments of our manuscript. Some of the questions raised were also noted by the other referees, which suggests that we did not fully explain certain key aspects as clearly as we should have. At the outset, it is important to highlight two key issues, as they underpin many of the subsequent questions and responses.
Firstly, we would like to emphasize that the aim of this study was not to conduct basin-scale hydrological balances or to assess runoff and flooding risks in specific areas. Rather, the primary purpose was to evaluate and better understand the role of several key parameters—commonly used in soil science as proxies for specific process related to water dynamic—in the assessment of soil-based hydrological ecosystem services and soil health, as illustrated in Figure 1. Accordingly, the study examines the impact on water regulation ecosystem services of parameters that are not all directly incorporated into the modeling, such as clay content and organic matter.
Secondly, by examining three markedly different pedo-climatic contexts, each characterized by substantial internal variability, our aim is to demonstrate that relying on single parameters to infer soil-based ecosystem services and soil health should be approached with caution. For example, a high Ksat value in the upper soil horizon (30 cm, as reported in the LUCAS database—the benchmark database for the European Commission) does not necessarily imply high infiltration. Infiltration is strongly influenced by the interaction among soil layering, boundary conditions, soil depth, and other factors. For instance, the presence of a compacted horizon below the topsoil, particularly in sandy soils with low water-holding capacity, combined with a skewed rainfall distribution, can significantly experience a not very-high infiltration despite high Ksat in the upper layer. To capture these interactions and provide robust estimates of water balance components under diverse conditions, we adopted a process based modelling approach based on the solution of the Richards’ equation applied to 315 soil profiles, representing approximately 50 soil types, across three case studies spanning a climatic gradient from central Europe to northern Tunisia.
Below we provide our detailed responses to your specific questions and comments.
Q: On the description of the modeling case: I somewhat lack a proper description of the modeling process, parameterization, etc.: a clear delineation and discussion of which of the reasonably important processes are (or are not) modelled by the model. Was the full functionality of the model used or did data limitation or some other reason reduce what processes were actually considered? For instance, preferential flow is a key process to consider given the aim of the paper. I also generally lack clarity on the parameterization of the model, how many profiles have been modeled, as well as any mention of validation of the model for the purpose or location(s).
A: For clarity and simplicity, we address the posed questions point by point:
Applied model: As correctly pointed out by Referee (and also by the Referee 2), the initial description of the FLOWS model lacked sufficient detail regarding its structure, underlying processes, and the relationships it employs. We appreciate this observation and have revised the manuscript to provide a more comprehensive explanation.
FLOWS is a physically-based hydrological model that simulates water movement within the soil profile. It is fundamentally grounded in the Richards’ equation, which governs unsaturated flow by accounting for the hydraulic properties of the soil and the gradients in water potential. This foundation ensures that the model captures the actual physical driving forces behind water redistribution, rather than relying on empirical or simplified assumptions. The model explicitly incorporates:
- Soil layering, allowing for heterogeneity in hydraulic properties across different horizons.
- Boundary conditions, both at the surface and at depth, which influence infiltration, percolation, and drainage.
- Flux-to-runoff processes, enabling the transition from subsurface flow to surface runoff under specific conditions.
- Dynamic redistribution of water, reflecting temporal changes in soil moisture due to precipitation, evapotranspiration, and internal flow
Since it is not feasible—also due to editorial policy—to provide a detailed explanation of all processes in a soil–plant–atmosphere water balance system, we will strengthen the description in the Materials and Methods section. In addition, we have published both the user manual and the theoretical background of the model, (10.5281/zenodo.17076522), which can provide the necessary details. Furthermore, As answered to the Referee 1, the model has been applied in several papers related to water and solute balance, even in complex soil profiles, both at the theoretical (see for example, Coppola et al., 2009; 2012; 2015a) and applicative levels (see for example Coppola et al., 2015b; 2019a, 2019b; Hassan et al., 2024; Wang et al., 2013; Porru et al., 2024). Therefore, we trust its robustness, and we are going to complete the references.
Given this physically based framework, strong criticism of the model's imperfections and the reliability of its simulations may not fully align with its intrinsic capabilities. While we recognize that the model does not capture all soil processes (e.g. lateral drainage), it is at the forefront of processing analysis (e.g. compared to the widely used bucket model). However, we acknowledge the limitations highlighted by the reviewer, particularly regarding the use of pedotransfer functions (PTFs). While PTFs are a practical tool for estimating soil hydraulic parameters from readily available data, they inherently introduce uncertainty due to their empirical nature and site-specific calibration requirements. Both the Reviewers and the Authors recognize these limitations, and we have taken care to discuss them transparently in the revised manuscript.
The issue of preferential flow is indeed complex and must be addressed in its various aspects. The model, currently in its beta version and not yet publicly deployed, incorporates preferential water and solute flow through a dual-porosity (dual-permeability) approach, enabling exchange processes between the matrix and macropore domains. In this configuration, the Durner model is used to describe both the water retention and hydraulic conductivity curves, reflecting a bimodal pore-size distribution.
However, at the scale of this study (simulation of 315 soil profiles in three areas across Europe and Africa), we do not have parameter values for each horizon of the individual profile, nor information on the presence or extent of macroporosity (even if limited to the surface horizon). In practice, while pedotransfer functions (PTFs)—despite their known limitations, as discussed later—are commonly applied and yield reasonable results under a matrix flow assumption using van Genuchten parameterization, no reliable PTFs currently exist for a dual-domain pore structure.
That said, we fully agree with the Reviewer that the issue of macroporosity and preferential flow is significant in understanding soil water dynamics. In response to this valuable observation, we have expanded our analysis by introducing an additional indicator to complement those already presented. Specifically, we now include a metric that estimates the nature and extent of macroporosity, calculated as the difference between saturated water content (θsat) and field capacity. This approach, inspired by the conceptual framework proposed in Nature [https://www.nature.com/articles/s41586-018-0463-x], provides a practical means of inferring the potential for preferential flow pathways within the soil matrix.
This addition enhances the comprehensiveness of our discussion and allows us to better address the complexity of water movement in structured soils. Accordingly, we will revise the Materials and Methods section to describe the new indicator and its rationale, and we will expand the Discussion to appropriately emphasise its implications for soil hydraulic behaviour and model interpretation.
Q: On the matter of scale: To me it feels that there is a lack of connection established between the scale of modeling and the scale at which conclusions are made: The model is apparently a 1D simulation model that has been run on individual soil profiles. However, “spatial differences” have been referred to (around L255), while there has not been any information provided on how the upscaling to 100s of km2 has been handled. Each area is obviously heterogeneous in terms of soil type and depth. I think this should be discussed, along with the limitations of the chosen approach.
A: If we have interpreted the question correctly, we agree with Reviewer 3 that no real upscaling was performed. Our intention, however, was to treat the profiles as independent from one another and as representative of a single soil map unit, thereby capturing all together the overall variability of soil properties within each study area. While this approach does not constitute a formal upscaling of processes, it does provide meaningful insight into the differences observed in each case study and, also, across the three. As suggested, in the limitations section, we will explicitly state that no real upscaling of the processes was carried out.
Moreover, the request for upscaling—together with other comments from Reviewers 2—convinces us that the main objective of the paper was not expressed as clearly as it should have been, which may have led to some misunderstanding. We will therefore revise the introduction and objectives to better highlight the central aim of the study and to avoid potential misinterpretation
Q: On using only the average year in the simulations: I find working with just a single “average” year to be a weakness of the study. Most of the events of concern – be that runoff, erosion, pollution, drought, etc. – are usually the products of either extreme events or extreme seasons/years. This is unfortunately not considered in the study, and I feel that the given advice is not based on the best available information. Working with a distribution of weather conditions (years) at each site would also allow the deduction of correlations between various factors in a much deeper way, allowing the study’s conclusions to be much better supported. Given that the data are available, the workflow would just need to be extended, and the results analysed at a higher dimensionality. I think the best value for the reader would be if the analysis was rerun on all available years and the correlation analysis was deeper, accounting for interactions and distributions that look into e.g. precipitation distribution (same soil under wetter and drier conditions) and soil depth as well, which is crucial when capacity-driven metrics are of concern.
A: We fully acknowledge the Reviewer’s observation regarding the nature of the processes investigated. As also noted in our response to Reviewer 2, our study was intentionally designed to explore soil behavior under average non-extreme conditions. This choice reflects the primary aim of the work: to evaluate and enhance our understanding of the influence of key soil parameters commonly employed in soil science to derive indicators of soil-based water-related ecosystem services (see Answer, p. 1), as illustrated in Figure 1. We recognize that this objective may not have been sufficiently emphasized in the original manuscript, and we will revise the text to make it clearer.
To be clear in the answer, the study does not aim to assess flood risk, nor does it address scenarios involving extreme events or climate change projections. These aspects, while highly relevant, fall outside the scope of the present work. Nonetheless, we appreciate the Reviewer’s suggestion regarding the potential for future investigations, particularly the use of hourly simulations to better capture short-term dynamics and extreme events. This is indeed a valuable direction for future research, and we will include a note to this effect in the revised manuscript.
Q: On the diversity of soil depths and their handling: Soil depth is obviously one of the keys to the findings, whereas it has not really been addressed as such. It has mentions, but it is a key factor that I believe requires much more attention. The study mentions e.g. drainage and ground water recharge. At what depth are they defined in this study? How are these factors handled when the soil is shallow?
A: According to the Reviewer's suggestions, we will expand the discussion on soil depth, which is one of the parameters considered in the assessment of ecosystem services and soil health. Regarding drainage and groundwater recharge, the two terms refer to the same process. Irrespective of soil thickness or groundwater table depth (if the groundwater does not intersect the soil profile), both describe the flux at the base of the profile—i.e., the outflow from the profile. In terms of ecosystem services, this flux corresponds to potential groundwater recharge. We recognize that this point may not have been sufficiently clear in the manuscript, and we will explain the concept in greater detail to avoid confusion.
Q: On some of the obvious uncertainties that are not being acknowledged and handled: Let me pick an example. I think the conclusions made regarding Ksat vs. infiltration at L417-429 are not thought through sufficiently and therefore are a bit careless. Essentially the issue is a set of limitations that originate from possibly several sources. In this study Ksat has been estimated by a generic PTF, which smoothens differences and is generally error prone. The other variable, infiltration is a result of modeling, which is prone to its own set of uncertainties and errors. Neither of them are measured here. Ksat and infiltration are not the culprits, but the lack of our ability to account for them is. Such sets of uncertainties should be acknowledged and discussed, and the findings placed in context, and the conclusions handled accordingly.
A: The Reviewer is certainly correct. While we did include these limitations in the original version of the manuscript, it is evident—also as noted by Reviewer 2—that they were not explored in sufficient depth. We will therefore expand this aspect more thoroughly in the revised version. In the spirit of the journal’s ‘open discussion,’ however, we wish to emphasize that the purpose and target audience of this work should not be overlooked. Our intention is not to move too far into detailed modeling analysis, but rather to remain anchored to the issues of soil health and the estimation of soil ecosystem services related to water regulation. We also consider it important to reflect on how the relevant scientific community currently approaches these topics.
As regards the Ks, however, would lead well beyond the scope of this paper: in addition to the uncertainties associated with pedotransfer functions (PTFs), there are also those linked to laboratory and field measurements (which can differ by one to three orders of magnitude for the same soil), the numerous simplifications required to approximate reliable pseudo-3D measurements in the field, and the strong variability of field data depending on the support size—even at the point scale. Therefore, as in many published studies, we accept that the use of PTFs—including the estimation of Ks—provides the best possible approximation at this scale and for the number of soil horizons considered (around 965 for all case studies). Our aim, therefore, was to highlight that the most widely available ‘dynamic’ hydrological measure in global databases must be used with caution. It is not a panacea that allows soil health to be estimated directly, without considering local processes and the evaluation of ecosystem services (as represented schematically in Figure 1).
Q: On advancing our understanding: The authors cite that past studies generally did not go to sufficient depth (L40-41), however, as a reader I do not really see how this study advances the state-of-the-art all that much. Much of the conclusions are rather common sense and reinforce relatively basic, textbook type knowledge that is relatable to soil depth and the “average” soil’s physical and hydrological characteristics. To refer to some examples, several very logical connections can be made in Table 2 already: infiltration is strongly influenced by Ksat (and likely soil thickness), runoff is strongly connectable with prevailing precipitation intensity and soil texture type, drainage by precipitation characteristics and Ksat. If there is no rain, there is no infiltration or runoff. If the soil is more course textured, there is good infiltration and no runoff. Examples are in e.g. lines 285-286, 291-318 and 336-337.
The quantitative confirmation by modeling is appreciated of course, but I feel that the reader does not receive sufficiently impactful new information that is of interest to THEM.
A: We appreciate the Reviewer’s comments and would like to clarify our position. While some of our findings confirm expectations that could be considered intuitive or self-evident, others appear to contradict what the reviewer refers to as 'common sense'. Furthermore, as mentioned previously, and as discussed extensively in the introduction, modeling should not be regarded as a purely didactic exercise but rather placed within the broader context of the soil science community and aligned with emerging frameworks such as the Soil Monitoring Law. These frameworks aim to assess soil health at the “soil district” level. It is precisely the prevailing assumptions regarding water regulation, ecosystem services, and the role of the soil–water nexus in soil health that we seek to challenge. Moreover — and this is far from a minor detail — the study evaluates the impact on ecosystem services of parameters not directly included in the modelling, such as clay content and organic matter.
A notable example is the relationship between crop water stress and organic matter content. It is commonly assumed that higher organic matter content increases water availability for crops, a notion supported by many studies. However, our results suggest that this does not necessarily result in a proportional reduction in crop water stress. As can be seen in Table 3 for the Marchfeld and Bologna sites, for example, the expected negative correlation is not observed. This discrepancy highlights the importance of distinguishing between water retention capacity and actual water storage, and therefore crop water stress, which is influenced by multiple interacting factors including rainfall distribution, soil structure, porosity, and layering.
Another example, in addition to the previously discussed Ks, is the case of AWC. Although AWC is widely used as the primary indicator of soil water storage capacity and crop water availability, it does not fully capture the complexity of the system, as our results demonstrate. Another relevant example concerns the notion that higher levels of organic matter enhance infiltration. However, as pointed out by the Reviewer (https://doi.org/10.2136/sssaj2004.0055), some studies have reported contradictory results — for example, a negative relationship between organic matter content and Ks in both newly developed and existing pedotransfer functions (PTFs).
Summing up, our results show that these relationships are dependent on soil and climate conditions. Therefore, they are not universally valid and must be interpreted within the broader context of system complexity. As requested, we will expand this discussion in the revised version.
Q: The authors should generally recognize and discuss the study’s limitations that may root from the model itself, the choices made in parameterization, the scale problems, simplifications (depth, climate “average”), and also data sources. I think the value of the reported study could be enhanced a lot by both addressing several of these with improved analysis, and discussing the rest.
Some detailed comments:
Q: L40-41: This statement invites that the authors make a very clear statement later on how their study improves on this.
A: We are going to change the Introduction as mentioned previously.
Q: L43-48: Here I recommend consulting e.g. https://doi.org/10.1016/j.earscirev.2021.103873 and https://www.nature.com/articles/s41586-018-0463-x and citing them if they are found relevant and fitting.
A: We are going to add the references to the Introduction as mentioned previously.
Q: L151-152: Please provide more on this – how the rule/decision was made which condition to apply. Also: was there any warm-up period considered for the model?
A: We are going to add all the requested details of model setting in the text.
Q: L155-156: How were those parameters obtained?
A: We focused on durum wheat because it is a widely cultivated rainfed crop in the three study areas, and its growth cycle is well-documented. The crop development stages were determined by integrating the FAO56 guidelines with expert knowledge, literature, and site-specific field observations. This approach ensured a robust and site-appropriate parameterization of crop development for each case study. For example, based on the typical crop calendar and the growth stages described in FAO56, flowering occurs earlier at RMEL than at Marchfeld and Bologna; while sowing takes place earlier at Marchfeld and Bologna compared to the Rmel site. To further strengthen transparency, we will provide the representative development stages for each site in the supplementary material and revise the main text accordingly.
Q: L167: How were the h1-h4 thresholds handled? Crop specific? Fixed? Please provide info.
A: h1-h4 thresholds were set as crop specific. We are going to add all the requested details in the text.
Q: L170: The choice of median year largely dismisses the significance of infiltration capacity in really wet years. Please acknowledge this appropriately.
A: As also pointed out by the Reviewers 1 and 2, we are going to acknowledge the limitations of our choice in the text. However, it must be stated that, in addition to selecting the median rainfall year for each site, we also examined the intra-annual rainfall distribution to ensure that the selected year was representative of the typical seasonal rainfall pattern of each site. Eventually, while the current analysis is based on median weather years to provide a baseline characterization, we are already working on extending the assessment to include both dry and wet extremes. This will allow us to better capture the variability and potential vulnerabilities of soil functions under climate variability and change.
Q: L185-186: Please use the original 1999 Geoderma reference for the HYPRES PTFs.
A: We are going to revise the reference according to the Reviewer’s suggestions.
Q: L207-210: At least some of the outputs and inputs/indicators are rather closely correlated. Could you please justify their handling? Example: infiltration vs. Ksat. Also: define AWC.
A: We are aware that some of the selected soil parameters are correlated with the model outputs. Their inclusion was intentional, as they are widely recognized physical and chemical indicators, as well as proxies, of water-related soil functions and processes. These indicators were limited to the upper soil horizon, in line with approaches adopted in several studies and projects (e.g., SERENA, BENCHMARKS) and consistent with the EU LUCAS database, which primarily targets the topsoil layer. The available water content (AWC) was applied according to its classical definition, as the difference between soil water content at field capacity and at the permanent wilting point (AWC = θFC – θWP), and this will be explicitly clarified in the revised text.
Q: L231-234: Please clarify the distinction between statistical and multivariate “relationships”. Latter is also a statistical method.
A: The Reviewer is right that multivariate relationships are also statistical in nature. In the revised version, we will clarify the terminology: by statistical relationships, we referred to simple bivariate associations (e.g., correlations or regressions between two variables), while multivariate relationships indicated more complex statistical approaches involving multiple predictors and responses (e.g., multivariate regression, PCA, or clustering). We will revise the sentences to make this distinction explicit and avoid ambiguity.
Q: L244: Is soil depth accounted for in the modeling approach? Please describe. Table 1 suggests that there have been huge differences in soil depth _within_ each study area. Soil depth of course, matters a lot when it comes to the soil’s buffering capacity in e.g., infiltration capacity, total pore volume in the profile, total available water for the crop, etc.
A: Yes, soil depth was accurately accounted for in all the 315 simulations carried out by the model. Moreover, we revised the Table 1 since some values were wrongly reported, even if input data were correctly considered in the simulations.
Q: L264-265: These differences in coarse and fine textured soil Ks values are barely different. Usually the variability of local measurements would not even make such values statistically different. I recommend double checking and also noting this in the manuscript. PTFs tend to smoothen within the available data range, and reduce such differences normally seen in the field.
A: Please, refer to the above answer A1 about the Ks and the PTFs.
Q: L272: “all” soil profiles? Please clarify what this really means. I also note that averaging over a huge diversity of soils within each region may simply mask the hot-spots, that matter a lot in reality.
A: “All soil profiles” refers to the total number of profiles simulated for each case study (Marchfeld: 215; Bologna: 73; RMEL: 27). Table 2 provides a summary of overall statistics to offer a comprehensive overview of the results. Nevertheless, we conducted a detailed examination of each individual profile to identify potential hotspots.
Q: L336-337: re: stress correlations: I think these were logically expected. Bologna has the deepest soils on average, as well as the most clayey ones, and the highest amount of precipitation. Increased clay just adds to retention, and it is likely that the “h1-h2” (wet) stress limit kicks in more often, so the reduction of clay has a positive effect. In Rmel the conditions are exactly the opposite: dry region, highest sand content, shallowest soils “on average”. Then the “h3-h4” (drought) stress limit is likely often met, so any clay present will help retain water and make it available to plants. To me this feels rather simple “on average”.
A: Yes, we fully agree with the Reviewer and will emphasize these aspects in the revised version of the work.
Q: L336-339: Why not run all years and then plot and correlate precipitation, clay content and CWSI? A cloud of points and deeper analysis of the potentially available data would lend much more support to the claims made.
A: The trade-off between model complexity, the number of soil profiles (n = 315) included in the study, and their discretization into 100 nodes for solving Richards’ equation necessitated restricting the investigation period to a single year, deemed representative of the climatic conditions of the three zones. In addition, the presentation of the results was challenging owing to the substantial number of parameters, processes, and cases considered.
Q: L340: This is against classic textbook knowledge, although it has been reported by others before. The authors could attempt to go deeper in explaining it. A starting point can be https://doi.org/10.2136/sssaj2004.0055 or https://doi.org/10.2136/vzj2016.03.0021, but there are a few other references available since.
A: We will address this and other related issues in greater detail in the review. As previously mentioned, the main objective of this work is to challenge the conventional assumptions often presented in textbooks and to account for system complexity when identifying relationships among parameters, indicators, and ecosystem services, with the ultimate goal of improving soil health assessment.
Q: L385-388: I think the main point is not to increase some correlation coefficients, but to help understand how variable interactions help explain why certain correlations between variables and indicators are seen and how they interact, especially in “non-average” cases. I suggest trying to focus on what can be generalized for the benefit of the independent reader, who do not associate directly with the cases of Bologna or Marchfield or Rmel.
A: The key point of this work is precisely not to generalize the results, but rather to carry out a site-specific analysis, depending on the combination of climate, soil, and—where relevant—management practices, which were not considered in this study.
Q: L392: What T values are referred to? Introduce the statistic used and show the data that you are referring to.
A: The t-values refer to the t-statistics used in hypothesis testing of each regression coefficient and the values of the t-statistic used to test the significance of each regression coefficient in the final stepwise regression model. The t-values in the regression output are the calculated t-statistics for each coefficient, obtained by dividing the estimated coefficient by its standard error. They are used to test whether each predictor (Clay, depth) has a statistically significant effect on the dependent variable (INF).
The t-values eported in the regression output (8 for Clay, and 2.8 for depth) are the values of the t-statistic used to test the significance of each regression coefficient in the final stepwise regression model, which indicating that clay has three-times level of confidence as a predictor.
We are going to add all the requested details to the text and change the capital T to t.
Q: L399-401: Without depth considered, basically the same thing has been stated already on L336-337.
A: We decided to discuss the contribution of each parameter to the correlation coefficient within the multivariate analysis, in order to improve readability. However, we will try to avoid unnecessary repetitions.
Q: L400-402: On L335-336 it has been stated that in Bologna moderate correlation was found and higher clay content has been associated with increased stress. Here it is explained that only depth (i.e. not clay) showed correlation with stress. So which one is it? Running different analyses should not be the explanation. Then one of them is the incorrect method.
A: The sentence “In contrast, at Bologna, this negative correlation is observed only for the soil depth variable” was related to the previous ones, where in the RMEL case study negative correlations were found with both soil depth and clay content. For the Bologna case study, instead, a moderate negative correlation was observed only with soil depth, while the correlation with clay content was positive. We will revise the sentence to improve readability.
Q: L416: The variability within each of the 3 areas has been greatly masked, whereas that is what would inform the reader how to judge their own cases. I don’t really see deep enough support for that in this study.
A: We are not fully sure to understand what is meant here by “The variability within each of the 3 areas has been greatly masked” as in each area soil variability (if this is what the Referee refers to), has been described based on the most complete and up to date soil information available in each area. Of course, there would be more to explore but in each area soil-landscape relationships have been deeply investigated and soil variability across landforms is known and fully accounted for the soil profiles considered for simulation. In terms of variability of soil responses to external drivers, in particular climatic drivers and to crop management drivers, we agree that referring to a median year and to a single crop has a masking effect on the full range of possible soil responses, but the goal of this simulation exercise was to compare average soil responses in terms of provision of regulating ES under a set of standard and comparable conditions and comparing these responses with a set of ES indicators.
Q: L417-429 (and also to some extent L430-439): Ksat is being singled out here, so I come to its defense to some extent. The authors should deepen their arguments at least where prior knowledge exists. Ksat of course is known to show great spatial variability that is very difficult to capture by PTFs that use only basic physical soil properties that have very limited relationship with soil structure. The estimates by PTFs are known to be smoothened (often a lot) compared to the original data support, and be often weak descriptors of the original data, no matter what is claimed about them. I ask the authors to look up the R2 value of e.g. the used HYPRES European PTFs – they will be surprised how low it is (around 0.2, as I remember). The mean Ksat values presented on L265-266 show such sign of smoothing. The problem is by far not only with Ksat and infiltration not being strong-enough correlated, but that in this study a PTF was used – and a rather generic one. Essentially, this study did not have Ksat to work with. It had some estimations of Ksat to work with – i.e. not necessarily the correct data. Second: I feel that the model’s imperfections, simplifications, simulated processes have not been thoroughly accounted for and discussed. The authors try to draw conclusions between an estimated variable and an indicator whose simulation is dependent on a multitude of imperfectly modeled processes, under “average” climatic conditions, and smeared over varied soil depths/types. I think all such uncertainties should be acknowledged, and at least some that are feasible to handle in a more detailed way should be improved. Relevant conclusions should of course be modified (softened?) accordingly.
A: The Reviewer is certainly correct in pointing out the limitations of relying on PTFs. This remains a fundamental challenge, especially when addressing problems at a regional rather than a local scale. As we answered in our response to Reviewer 1 and 2, we will strengthen the discussion section to more explicitly acknowledge and elaborate on the limitations associated with the use of PTFs, so that the reader is fully aware of their implications in the present context.
As for the second part of the Reviewer’s comment, the Authors respectfully disagree with it, as the lack of a clear description of the model in the manuscript may have contributed to a misunderstanding. FLOWS is a physically-based model, grounded in the Richards’ equation, which explicitly accounts for the actual driving forces, soil layering, boundary conditions, flux-to-runoff processes, and the dynamic redistribution of water within the soil profile. The strong statements concerning the model’s imperfections and simulations do not align with the intrinsic nature of the model, while we nevertheless acknowledge the limitations highlighted by the Reviewer, as well as by the Authors themselves, regarding the use of PTFs.
Q: L440-455: I feel that the authors could have done a deeper analysis to support this section with more results, and make it more concrete. Rainfall distribution could have been characterized, the frequency of intense (defined as…?) rain events quantified, etc.
A: The Authors will provide additional details on the distribution and characteristics of rainfall to further clarify this discussion.
Q: 456: This is exactly what is needed – basing the study on appropriate data. I don’t want to say that a study based on PTF estimates is illegitimate, but the conclusions drawn should account for the source of errors that originate from PTF use, model use, spatial and temporal simplifications and averaging, etc.
A: We are going to highlight the study limitations as requested by the Reviewer.
Citation: https://doi.org/10.5194/egusphere-2025-1927-AC3
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 715 | 72 | 34 | 821 | 18 | 30 |
- HTML: 715
- PDF: 72
- XML: 34
- Total: 821
- BibTeX: 18
- EndNote: 30
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
General comments
The authors studied the link between soil physical and chemical parameters and indicators of hydrological processes and ecosystem services. They performed an extensive simulation study with the hydrological FLOWS model for 3 contrasting regions regarding soil variability and climate. They showed that the most relevant soil properties for soil functions and ecosystem services depend on the environmental conditions in a region, which is an important innovation.
The paper reads fluently, the methodological setup is sound, and the results are well presented in Tables. However, I still have questions and remarks, as described below.
Particular comments
Figure 2: For each site, a different soil legend is used. This makes comparison hard. Can you create a figure with a single uniform soil legend?
Figure 3: In the Figure of FLOWS, you show the groundwater. Does FLOWS consider capillary rise?
Also related to this Figure: “lateral drainage” is drawn in the unsaturated zone. As the majority of lateral drainage occurs in the saturated zone, better to depict “lateral drainage” there.
Line 130: The FLOWS model plays an important role in your analysis. Can you mention studies that showed the reliability of FLOWS?
Line 147: At which soil depth did you apply the free drainage option?
Line 159: You define here the potential soil evaporation, which in general is much larger than the actual soil evaporation, an important part of your water balance. How did you calculate the actual soil evaporation?
Line 165: Is potential root uptake also assumed to be proportional to root density?
Lines 174-183: You select for each site a median year based on the amount of rainfall. Did you also check whether the rainfall distribution over the year was representative of the site's climate?
Precipitation input: At which interval did you provide input of rainfall fluxes? In order to simulate surface runoff, daily rainfall intensities are generally not suitable.
Lines 184 and 185: The data of soil texture and organic matter: to which soil layer do they belong? Did you assume these data to be uniform over the entire root zone depth (max 80 cm)?
Line 187: The crop coefficient is defined depending on the stage of crop development. How did you derive these stages?
Table 2: The water balances of Bologna and Rmel do more or less close. However, the fields in Marchfeld, on average gained 94 mm of water. This is a substantial amount compared to other water balance components and a bit strange for a median year. Please clarify this in the text.
Line 274: Minimum and maximum values are not depicted in Table 2.
Lines 417 – 433: Important that you mention the ambiguity on the role of the saturated conductivity. You might add here that the simulated amount of runoff not only depends on the saturated conductivity but also on the time intervals with which rainfall is provided. In general, daily intervals underestimate the runoff amounts.
4 Discussion: Mention important limitations of your study:
Detailed comments
Lines 110 – 114: Current text is wrongly phrased and repetitive. Define the climate of site Rmel more concisely and clearly.
Line 134: “Genucthen” should be “Genuchten”.
Line 138: The symbol “h” was already defined in line 134. Define each variable once, upon first usage in the text.
Lines 189 and 190: Include references.
Line 191: Include units.
Line 134: “groupf” should be “group”.
Table 2: In the caption, clarify CWSI.
Table 3: In the caption, clarify SD.