Effects of passive storage on modelling hydrological function and isotope dynamics in a karst flow system in southwest China

Abstract. Representing passive storage in coupled flow-isotope models can facilitate simulation of mixing and retardation effects on tracer transport in many natural systems, such as catchments or rivers. However, the effectiveness of incorporating passive storages in models of complex karst flow systems remains poorly understood. In this study, we developed a coupled flow-isotope model that conceptually represents both “fast” and “slow” flow processes in heterogeneous aquifers to represent hydrological connections between hillslopes and low-lying depression units in cockpit karst landscapes. As this model originally included a varying number of passive storages at different positions of the flow system (e.g. fast/slow flow reservoirs combined with different hillslope/depression units), the model structure and relevant parameters were optimized using a multi-objective optimization algorithm. This was used to match detailed observational data of hydrological processes and isotope concentration in the Chenqi catchment in southwest China. Results show that the optimal structure for a coupled flow-isotope model incorporated only two passive storages in fast flow and slow flow paths of the hillslope unit. Using fewer or greater numbers of passive stores in the model could lead to under- or over-mixing of isotope signatures. This optimized model structure could effectively improve simulation accuracies for outlet discharge and isotope signatures, with > 0.65 of the modified Kling-Gupta efficiency. Additionally, the optimal tracer-aided model yields reasonable parameter values and estimations of hydrological components (e.g. more than 80 % of fast flow in the total discharge). Furthermore, results imply that the solute transport is primarily controlled by advection and hydrodynamic dispersion in steep hillslope unit, which is a remarkable phenomenon in the karst flow system. The study resulted in new insights, more realistic catchment conceptualizations and improved model formulation.


also exists in Southeast Asia, Central America and the Caribbean. Our selected catchment of Chenqi is a karst experimental catchment focused on investigations of hydrological, ecological and geological (carbonate dissolution) changes under climate change and human activities. So there are detailed observational data and field investigations in this catchment. The relevant publications cited in this study are necessary to provide the background context of our new model development and analysis. Although Chenqi catchment is small, the geomorphologic characteristics can represent a broad region of headwater catchments in cockpit karst landscapes.
In the polje/tower karst areas, the depression is more interconnected with isolated towers scattered throughout the terrain (Lyew et al., 2009). Geological surveys and observations show the hillslope unit lacks surface flow, and the depression unit has surface and underground drainage networks in such karst areas, including our study catchment. Understanding of interconnections of flow systems are vital for developing conceptual hydrological models for cockpit karst landscapes. Our new model presented in the paper is based on the coupled hydrology-isotope model developed by Zhang et al. (2019), a co-author of this manuscript. In this earlier model, the cockpit karst catchment was divided into two morphological units (hillslope and depression) and three water storage compartments (reservoirs) (hillslope reservoir, fast flow and slow flow reservoirs in depression). We substantially improved the model structure with a binary flow system (fast flow and slow flow ) in the hillslope unit, and the functioning of a binary moisture storage system of unsaturated zone (see Fig 4 in the original manuscript). Moreover, we optimized the model structure with a varying number of passive storages at different positions of the flow system (e.g. fast/slow flow reservoirs combined with different hillslope/depression units) based on a multi-objective optimization algorithm for best matching detailed observational data of hydrological processes and isotope concentration in the Chenqi catchment.
We agree there are various connections between hillslope and depression fast/slow flow reservoirs, and the model structure can be further improved in terms of the geomorphological surveys of the catchment. So, we set another reasonable connection between hillslope and depression fast flow and slow flow systems, and re-calibrated and validated the model (see descriptions below). We referenced the previous investigated results and will show more detailed geomorphological data in the revision (e.g. electrical resistivity tomography (ERT) image in Fig S1) to show how data has informed the evolution of this new model. Figure S1. ERT image in the study depression. They interpret the ERT results as (a) an upper layer consisting of moist soils or extensively fractured rock (marked in blue); (b) carbonate rock with a high secondary porosity (and hence permeability; marked in light blue/yellow); (c) an underlying carbonate rock with low secondary porosity and hence relatively low permeability (marked in red) (Chen et al., 2018).
Additional reference: Lyew-Ayee, P., Viles, H, A., Tucker, G, E.: The use of GIS-based digital morphometric techniques in the study of cockpit karst, Earth Surf. Process. Landforms., 32, 165-179, https://doi.org/10.1002/esp.1399, 2009 (1) Lines 60-65: the list of references could be extended by some more recent articles (<10 years) Reply: We carefully read the two references suggested by the reviewer, and added associated publications in the most recent 10 years as follows: The residence time: Brki, Z., Kuhta, M., Hunjak T.: Groundwater flow mechanism in the well-developed karst aquifer system in the western Croatia: Insights from spring discharge and water isotopes, CATENA., 161, (2) Lines 110-111: the question here is answered in the specific case of the study site. In what way is the structure of the model finally proposed transposable elsewhere and at a different scale? Are all cockpit systems of the same nature in terms of their hydrological functioning? Reply: Lines 110-111 "Particularly, the effects of passive storage structures are underexplored in terms of the location and number of passive storages needed for fast and/or slow flow reservoirs in hillslope and/or depression units, respectively. Consequently, it remains unclear what is the most efficient way of incorporating passive storage into coupled flow-tracer simulations." In the original manuscript, we focused optimization of passive storages in hillslope and depression fast/flow reservoirs. We set fourteen schemes (scenarios) that incorporate 0~4 passive storages into different positions within the karst flow system, i.e., fast and/or slow flow reservoirs in combination with the hillslope and/or depression units (Table 3). We obtained the optimal structure (model f) for the coupled flowisotope model that incorporated two passive storages in fast flow and slow flow paths of the hillslope unit. This optimal structure is obtained based on the hydrological connections of hillslope -depression fast flow (HF-DF) and hillslope -depression slow flow (HS-DS). We further consider another possible connection of hillslope -depression fast flow (HF-DF), and hillslope slow flow (HS)-depression fast/slow flow (HS-DF/DS) with a ratio of rhd of HS contributing to DS (Fig. S2). The optimized rhd is 0.39. It means that about 61% of hillslope slow flow can enter depression fast flow reservoir. The optimal model structure of the passive-active storage connections is the same as the previous result (model f in Fig. S2) while the optimized parameter values and hydrological components have some differences (See Table S2~S5). From a geomorphological aspect, in the polje/tower karst areas, the depression is more interconnected with isolated towers scattered throughout the terrain (Trudgill, 1985). We proposed the concept of "hillslope-depression-stream" continuum that can capture the morphologic features of the cockpit systems (Chen et al., 2018;Zhang et al., 2020a). So, the developed model was based on the spatial discretization of hillslope and depression units, each with characteristics dominating runoff generation processes, streamflow processes and hydrological connectivity. In our study, the runoff generation is estimated based on water balance in unsaturated zone storage, the streamflow processes are routed by hillslope and depression fast and slow flow reservoirs, and hydrological connectivity includes connections of unsaturated zone (recharge)-rhd saturated zone (storage), and hillslope (fast and slow) flow -depression (fast and slow) flow. We believe that this model structure captures the internal catchment processes and hydrologic pathways of cockpit systems. Since hillslope runoff is regarded as a "water tower" for supplying the depression agriculture. Understanding the hillslope and depression hydrological functionality and their connections is necessary.
In our model, we used a distribution curve of the unsaturated storage capacity to describe the spatial heterogeneity of storage volumes, and fast flow and slow flow systems to elucidate dual karst flow system on a large scale (e.g. hillslope and depression units). Such delineations have been proven to be effective in other conceptual models, such as the VarKarst model (Hartmann et al., 2013) and the Xinanjiang model (Zhao, 1992). Surely, the model parameters still need to be calibrated when the model is applied to other catchments, but in principle the modeling approach is transferable. In large catchments, the model should incorporate river and channel routings that can play an important role in streamflow variations. (3) Study area section: The description of the site, especially the depression area, is very small. As expected, the average soil thickness is lower in the hillslope unit than in the depression unit. It is also expected that the nature of the soils is different and therefore also the field capacity. Can the authors provide details on these field characteristics? Also, please explain the phrase "perennially flowing underground conduit connecting the hillslopes to the catchment outlet" (see also line 337-338). Do you mean that there is a main karst conduit within the depression that transmits water to the outlet? Reply: We have done lots of in situ sampling and analysis as well as field surveys by using electrical resistivity tomography (ERT) in the study catchment, particularly in the flat depression (see Fig S1). In depression, the accumulated soils are thick (~200 cm) and cultivated for crops of corn and rice paddy. The soils are most silt loam consisting of over 80% of clay and silt with soil particle size of smaller than 0.02mm and bulk density of 1.31 g/cm 3 . The soil porosity ranges from 32% to 47%. In the hillslopes, Quaternary soils are thin (less than 30 cm) and irregularly developed on carbonate rocks. Outcrops of carbonate rocks cover 10%~30% of the hillslope area. The soil at sites from the shallow to deep layers in the catchment varies from sand loam, consisting of mostly sand (56~80%) and fine sand (20~40%), to calcareous soil and silt (1~10%). The bulk densities increased with depth, ranging from 1.02 to 1.33 g/cm 3 (Chen et al., 2009).
The depression unit has an underground channel/conduit system with perennial flow (see blue color in Fig S1). The high permeability zones (conduits) are sporadically distributed at the upper depression (the hillslope foot), which collects the hillslope flow. These connections can be identified by the slow recession of water table for the well (W4) at the foot of hillslope after rainfall ceases (Chen et al., 2018). The widely distributed conduits in the upper depression are gradually concentrated to an underground channel at the catchment outlet.

(4) What is there between the 2 m depth at the base of the soil and the water level reached at 13-30 m (see line 187)? How deep is the bedrock? What is its nature?
What is the nature of the water table in the depression? It is doubtful that we are still in a karst system of the same nature as the hillslope.

Reply:
The depression aquifer consists of the soil layer (about 2m thickness) overlying the lower fractured rocks according to the ERT image (see blue and light blue/yellow in Figure S1) and the drill core sampling (see Fig 3 in Chen et al., 2018). So, the depression aquifer has the bedrock (the impervious marlite formation) at depths of 30~50 m.
On line 187, the depths of 13 ~ 35 m below ground surface refer to the sampling depth of groundwater, instead of the water levels.
The depression is located in the low-elevation area (<1340 m) and steeper hillslopes have high elevations ranging from 1340 to 1500 m as shown in Table 1. The water level ranges are 1,267.4~1,275.9 m at W1 with a mean of 1273 m, and 1,280.0~1,285.2 m at W4 with a mean of 1282 m. There is a main underground channel in depression with an ascending spring at the catchment outlet, and high flows can spill over bottom of the depression ditches (referring to surface river channel with overland surface flow in Fig S2). So, in Fig 1, the two points at the outlet refer to the observation sites of underground and surface river channels at the catchment outlet (see Fig S3). The discharge used for simulations is the total of subsurface and surface discharge (see Fig S4). The discharge of surface flow and underground flow is 45% and 55% of the total discharge, respectively.
Two hillslope springs can be observed in the study catchment. We selected a perennial spring at the hillslope foot. The location has been added in the figure (see Fig  S3).  (6) Lines 164-166: Do the hydrographs mentioned refer to those observed at the outlet? There is ambiguity because the following sentence refers to epikarst springs.

Reply:
Lines 164-166: The hydrographs refer to the total discharge of surface and subsurface streamflow (see Fig S4).
The hillslope springs are formed by an impermeable layer (marlite) underlying the fracture zone (epikarst). We will change epikarst springs to hillslope spring to avoid misunderstanding.
Observational dataset: (7) The interpretation of figure 2 is very questionable, especially the regression line of the W1 points. There is a very large dispersion but the points W1, W4 and hillslope spring are to be included in the same O/D relationship which is also that of the local rainfall. The purple, green and black lines are disturbed by the few points indicating evaporation. The grey points (outlet) under the rainfall line are divided into two groups, probably indicating an evaporation process under different relative humidity conditions. In my opinion, it is not possible to argue about the age of the water from isotopic enrichment or depletion information alone: 1) The differences between the means for each set are modest and the number of measurement points is different each time. The difference in the mean between W1 and W4 is of the order of the measurement error. These differences are therefore not statistically interpretable. 2) Can't the apparent enrichment of W1 and outlet (vs W4) come from the inclusion of evaporated water? the apparent enrichment of W1 and outlet 3) Why is the dispersion on W1 the lowest? Could there be a different origin of water in W1 and W4? (linked to the organisation of the fracture network in hillslope). 4) Also specify the number of points and the origin of the data for the LWML.

Reply:
We agree that the enriched groundwater is caused by evaporative isotopic fractionation. After checking the data points, we find that the two groups of the more enrichments of δD and δ 18 O with δ 18 O>-7‰ in Fig S5 occur in different periods for the outlet streamflow and hillslope spring. The upper points occur in the period of late spring ~early summer (May~ early June) with consecutive occurrence of small rainfall events while the lower points mostly occur in some summer days (middle July and August) after a long period without rain. In the period of late spring ~early summer, rain water is more enriched (see Fig. 7 in the original manuscript), resulting from the Westerly water vapor with low humidity and local moisture recycling, and thus the enriched rainfall infiltration controls isotopic concentrations of the hillslope spring and the outlet streamflow. By contrast, in the large rainfall period of summer, rainwater is depletion (see Fig. 7 in the original manuscript) due to the Western Pacific water vapor source with high humidity. Nevertheless, evaporation is strong in the dry period of summer. As streamflow recesses rapidly in the catchment, the low flow after a long drought period mostly comes from aquifer storage with strong evaporative isotopic fractionation.
The differences in δD and δ 18 O values at sites are related to the extent of mixture with new water (e.g., from rainfall recharge). For W1 close to the catchment outlet, the site is located at a locally confined aquifer underlain by rocks with poor permeability according to the ERT survey. So, the subsurface flow seldom mixes with new water (rainfall) (Chen et al., 2018), resulting in more enrichment of groundwater at W1.
The plot of the δD and δ 18 O relationship is not directly related to water age. We will delete this description. Figure S5. Plot of  18 O-D for catchment outlet discharge and hillslope spring (1) In Table 2, the sampling time is the same (wet season) for W1 and W4 although the sampling frequency is less than that of the hillslope spring and outlet discharge. The differences of the mean δD and δ 18 O values between W1 and W4 are larger than the measurement error (± 0.5 ‰ for δD and ± 0.1 ‰ for δ 18 O). To be comparable, we recalculated the statistical values using the data in the same period (Table S1). Although the mean and range of δD and δ 18 O values in Table S1 are different to those in the previous calculation in Table 2, they lead to the same conclusions as the previous result. (2) The depression groundwater at W1 and outlet flow is more enriched compared to that at W4, attributable to both evaporation and mixing with the new water (e.g. rainfall recharge). As W4 is located at the hillslope foot, and groundwater there receives more new water (fast flow) from hillslope; Table 2 shows the mean δD and δ 18 O values of W4 are closer to those of the hillslope spring and rainfall. W1 is located at a locally confined aquifer surrounded by rocks with poor permeability, and the flow seldom mixes with new water (rainfall) (Chen et al., 2018).
On the other hand, depression groundwater partly comes from rainfall infiltration and percolation through the thick soils (Zhang et al., 2019), which undergoes evaporative fractionation. Our re-optimized coefficient of the evaporative fractionation also supports this conclusion. As shown in Table S3 (model f), the coefficient of evaporative fractionation in depression (lsd = 0.05) is greater than that of hillslope (lsh = 0.01).
(3) The smallest coefficient of variation (CVs) is resulted from the less seasonal fluctuation of water table at W1 due to little rainfall recharge from the upper low permeability layer (refer to Figs 4 and 5 of Chen et al. (2018)).
(4) The LMWL of δD =8.18δ 18 O+9.52 comes from the daily rainfall sampled over the whole study period at Chenqi catchment. The number of points (253) for the LWML data is listed in Table 2. We have added a plot of the δD ~ δ 18 O data and the fitted line in Figure 2 (see Fig S6). Figure S6. Plot of  18 O-D for rainwater, catchment outlet discharge, hillslope spring and depression groundwater at wells W1 and W4 Conceptual model structure : (8) In connection with the study area section, one level of explanation is again missing for a satisfactory understanding of the system. The structure of the model logically foresees a dual flow system in the 2 units and in the 2 compartments ZNS-ZS. But it seems to me that the authors make two very strong assumptions that need to be justified: 1) The fast and slow reservoirs are in perfect connection between the unsaturated zone and the saturated zone. This can be understood in the karstic part of the hillslope system but the continuity does not seem so obvious in the depression part where the nature of the slow flow/fast flow partition can be quite different between the soil and the water table.
Reply: As shown in Fig 4 (and Fig S2), over most of the catchment area (i.e., the low permeability area of ), runoff generated (free water R) in the unsaturated zone connects with both slow flow (ks) and fast flow (1-ks) reservoirs of the saturated zone, in which ks is a discount coefficient of R entry into the slow reservoir ( Fig S2). In the remaining area comprising the high permeability area (1-), rainfall directly enters underground channels through surface-connected sinkholes commonly found in carbonate aquifers.
In the depression unit, there are still some sinkholes that can accommodate rainwater even though the coverage ratio is small (e.g., 1% according to the re-optimized parameter of (1-αd) for model f in Table S3).
(9) 2) The authors suppose a hydrological continuity between the slow and fast flowing reservoirs of the 2 units (see also lines 330-331). Are there any tangible arguments to assume that slow flows from hillslope will retain this slow flow property in the depression (same for fast flows)? Reply: Various lines of evidence have demonstrated the hillslope-depression fast flow connection. For heavy rainfall events, the observed hydrographs are primarily dominated by fast flow. In the mid-season after an extremely heavy rainfall, hillslope flow is highly synchronized with outlet flow, and the relationship between hillslope spring discharge and outlet discharge approaches a monotonic function (R.R. Zhang et al., 2020a). Figure S5 also shows that data of  18 O-D for the outlet discharge are strongly overlapping with those of hillslope spring for large rainfall events.
We agree that the hydrological connections of hillslope slow flow (HS) -depression fast/slow flow (DF/DS) are not perfectly conceptualised in our previous model structure (see Fig 4 in the original manuscript). When the depression water level is low and the storage deficit is high, a portion of HS could be concentrated into the depression conduits (DF). So, we redesigned connections of the flow system of the two units as shown in Fig S2 (i.e., HF -DF and HS-DF/DS connections). In this flow system, a parameter rhd is used to represent allocations of HS between DF and DS (Fig S2 and  Table S2). The optimized value of rhd is 0.39 for model f (Table S3). Other associated variables are correspondingly recalculated, as shown in Table S4.    In this context, the sentence referring to fast flows speaks of large fractures vs. swallow holes, which suggests that the latter formations are in the depression part. Can you confirm this impression? If so, should this be linked to the "perennially flowing underground conduit" mentioned in the study area section and the "underground channel in depression" in line 337-338? Overall, the authors should make an effort to describe the hydrogeomorphological context of the system and better relate this information to the structure of the proposed model.

Reply:
Your descriptions are correct. We will revise this part of the description and strengthen the geomorphological conceptualisation in the study catchment as described above.

(11) 3.1.2 isotopic concentration routing: you do not take into account isotopic fractionation, whereas Figure 2 shows that there is evaporation. This fractionation is however integrated in the model proposed in Zhang et al (2019). Can you explain why you chose to ignore this process? Reply:
Based on your suggestion and the model proposed by Zhang et al. (2019), our we will developed out model further and include isotopic fractionation by adding an isotopic fractionation coefficient ls. Then the mass balance in the unsaturated zone storage can be expressed as: where WU (WU=W+Wpas) is the moisture storage consisting of active storage W and passive storage Wpas, δp and δb are the stable isotope concentrations of rainwater (P) and moisture (and water yield R), respectively, and is the coefficient of evaporative fractionation.
As expected, the optimized ls value is larger in the depression unit (lsd =0.05) than in the hillslope unit (lsh =0.01) (see model f for Table S3). The parameter ranges are set relying on the physical characteristics established in previous studies in this catchment. Some parameter values are directly specified according to field investigations (e.g., the ratio of matrix flow area α) and suggestions by other studies (e.g., the storage capacity of moisture wm and the passive storage for slow flow and fast flow Vpas) (Xue et al., 2019;He et al., 2019;Liu et al., 2020).
We did further calibration for all parameters. The recalibrated parameter values are different to those of previous calibrations, but the parameter values in the hillslope and depression units are ordered similarly to the previous results. For example, the ratio of matrix flow area of α in hillslope is 0.94, smaller than that in depression (0.99); wm representing the soil moisture retention capacity of 56 mm for thin soils over hillslope, is much smaller than 82 mm for thick soils over depression (see the parameter values for model f in Table S3). b represents spatial heterogeneity of water storage capacity for the matrix of unsaturated zone, instead of the conduits (fast flow) because the conduit area has been separated from the area for each of the hillslope and depression units (see Fig S2).
The curves of storage capacity WM related to the proportion (f/F) of the matrix area for unsaturated zone are shown in Fig S7. It shows over half of the area with WM less than 62.2mm and 94.3mm for thin soil in the hillslope unit and thick soil in the depression unit, respectively.

Figure S7. Storage capacity curve
Most studies show that the volumes of passive storages (Wh,pas, Wd,pas,Vs,pas and Vf,pas in Table S2) are generally one order of magnitude larger than those of active storage (Dunn et al., 2010, Soulsby et al., 2011, Ala-Aho et al, 2017. When the parameters are calibrated in this study, the calibrated values are 535 and 509mm for the unsaturated zone passive storages (Wh,pas and Wd,pas) in the hillslope and depression units, respectively, and 316 and 325mm for the saturated zone passive storages (Vfh,pas and Vfd,pas) in the hillslope and depression units, respectively (see model f in Table S3). These large passive storages ensure damping and time-lags of  18 O and D in streamflow response compared with precipitation fluxes, implying large mixing volumes that are usually much greater than dynamic storage changes estimated by water balance calculations (Birkel et al., 2011;Fenicia et al., 2010;Soulsby et al., 2011).
(13) In the original manuscript, model n includes all the parameters associated with passive storage for calibration. Other models include different choices of passive storages for fast and slow flow reservoirs in the hillslope and depression units. For example, model f includes two passive storages of hillslope fast and slow reservoirs as the optimal model structure excludes other two passive storages of depression fast and slow reservoirs. So, the parameters of exchange coefficients between active and passive storage for the depression unit (sd and fd for model f in Table 4) are not necessary.
"NA" in Table 4 means not applicable. We will use "-" to replace "NA" in order to avoid misunderstanding.
(14) Line 383: prefer μ to σ to express an average (or a ratio of averages) Reply: We will change the expressions of the two variables.
(15) Line 397: I am not very familiar with multi-objective optimisation algorithms but the number of iterations seems low. In fenicia et al (2007), the number of iterations is rather in the order of a few thousand.

Reply:
We have executed the multi-objective optimisation algorithm again by increasing the number of iterations from 100 to 1000. The calibration of the NSGA-II algorithm was performed as follows: For a number of iterations (e.g. 1000 in this study), 50 parameter sets were initially retained for screening out the sets with Biasq less than or equal to 0.2. Then the remaining sets are sorted from the largest to the smallest according to the sum of corresponding KGEq and KGEc. Finally, 30 sets are selected as the Pareto-optimal solution (Nan et al., 2021). The corresponding objective function values (average of the optimal solution sets) for both the calibration and validation periods were extracted.
However, even with the number of iterations increased from 100 to 1000, we find that the results are not significantly different.

Performance of models:
(16) Lines 428-429 (and Figures 6-7): The authors should be more critical about their results. In particular, the model does not really succeed in capturing isotopic variations. In many cases, it overestimates or underestimates the observed values, especially in the calibration period. In validation, the results are better because there is less variability. Finally, the performances are not better (or even worse) than those obtained with the model of Zhang et al. How can these shortcomings be explained in terms of the structure (defaults) of the model? The announced contributions of the fast flow reservoir are not inconsistent, but it is not reasonable to justify these results by those of another model. Once again, there is a lack of arguments from the field.

Reply:
For the coupled modeling of hydrological and chemical or isotopic processes, all models give a higher accuracy of streamflow simulation than chemical or isotopic simulations (Soulsby et al., 2015;Dehaspe et al., 2018;Mudarra et al., 2019;Birkel et al., 2020). Additionally, most isotopic simulations of the coupled models are compared with daily or monthly analysis data of the isotopic concentrations, and seldom compared with exacting test of hourly data as executed in our study. Our simulated accuracy is acceptable for the isotopic process with the KGEc of 0.59 and 0.73 in the calibration and validation periods, respectively, compared with that from other studies (Delavau et al., 2017;Aaron et al., 2019). The KGEc is also greater than 0.5 from the simulation of Zhang et al. (2019) in wet season.
The optimized model (model f) captures the sharp rise and decline of high flow and isotopic variations, but it can not simulate some low flow and isotopic processes similarly well. Particularly, the model overestimates some fast decreases of low flow caused by groundwater pumping for agriculture use (i.e., the sudden declines of streamflow and isotopic concentrations in June as shown in Figs S8 and S9). Comparatively, the validation period has fewer low flows and less groundwater pumping influences. As a result, the simulations are better in the calibration period than in the validation period.
We agree there are some limitations in our developed conceptual model due to the strong heterogeneity of the karst media. For example, some runoff can occur from areas of impermeable rocks despite small rainfall. We will discuss the limitations of the model in more detail in the revised manuscript.
We will strengthen descriptions of the field investigations used for comparison with our simulated results. As listed in Table S4, the simulation results of model f suggested that the proportion of the total subsurface flow (slow flow and fast flow at underground channel) is 58%, and surface flow from the surface channel is 42% of the total catchment flow. These proportions are consistent with 55% and 45% from the observations at underground conduit outlet and surface channel outlet, respectively.  Tracer-aided modeling in the low-relief, wet-dry tropics suggests water ages and (17) Table 5: How can it be explained that the parameters for the validation are better than those for the calibration? Did the authors consider switching the 2 periods? Reply: Please refer to the above explanations why the simulated discharge and isotopic variations in the validation period are better than those in the calibration period. Since the calibration period is long and includes different magnitudes and variations of hydrographs, the calibrated parameters using the data are more representative than those from the data in the calibration period.
On the other hand, not all models obtain a higher accuracy of simulations in the calibration period. For example, models d and h in Table S5 obtain a lower simulation accuracy in the calibration period, compared to that in the validation period. These two models neglect passive storage of hillslope fast flow reservoir while they include passive storage of hillslope/depression slow flow reservoirs. No passive storage in the hillslope fast flow reservoir reflected the input signal overestimating stream isotopes. This was also reported by Birkel et al. (2011) for mountainous catchments in Scotland. We will concisely describe this as follows: H and D represent hillslope and depression units, respectively; F and S represent fast flow and slow flow, respectively; P represents passive storage. Thus, the dual flow system (F and S) combining with the two units (H and D) is represented by HF and HS, and DF and DS, indicating hillslope fast flow and slow flow, and depression fast flow and slow flow, respectively.
The flow system consists of hillslope and depression fast flow connection (i.e.. HF-DF), and hillslope slow flow (HS) and depression fast/slow flow connections (HS-DS/DF).
The passive storage (P) added in each flow reservoir (HF, DF, HS and DS) is represented by the subscript P, i.e., HFP, DFP, HSP and DSP.
So, model a ~n could be described by: