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 

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). 

1. 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.  

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 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).