Enhancing process interpretation with isotopes: potential discharge-isotope trade-offs in ecohydrological modelling of heavily managed lowland catchments

Tracer-aided modelling (TAM) can help quantify ecohydrological processes, as stable water isotopes (δ¹⁸O and δ²H) can provide complementary information beyond streamflow and help constrain equifinality. Whilst TAM has been successfully undertaken in smaller rural research (<100 km²) catchments with limited anthropogenic impacts, its utility in more heavily managed catchments remains untested, particularly as isotope samples are usually unavailable. This study investigated four sub-catchments (Berste, Wudritz, Vetschauer, and Dobra) in the heavily-managed Middle Spree River basin (ca. 2800 km²), in NE Germany; a strategically vital water resource supplying drinking water to the capital of Berlin, Germany, and sustaining agricultural and industrial demands. Disentangling ecohydrological water partitioning in this evapotranspiration (ET) dominated region is complicated by heterogeneous land use, extensive hydraulic infrastructure and a long legacy of intensive management. We used the spatially distributed tracer-aided model STARR to simulate water fluxes and storage dynamics over a 6-year period. Temporally coarse (seasonal) isotope data was used in calibration as well as streamflow to help constrain estimates of subsurface runoff sources and ET partitioning. In most cases balanced calibrations using both isotopes and discharge increased confidence in plausible process representation in the modelling. When the trade-offs between dual calibration targets could not be reconciled, these were likely explained by anthropogenic factors that were not easily incorporated in the modelling framework. Such trade-offs therefore provide opportunities for learning about epistemic errors (e.g. un-represented water withdrawls for irrigation) that can be used to improve future models for heavily managed catchments. Our modelling framework shows the potential for informative insights from wider use of even sparse isotope data sets in tracer-aided modelling of complex, heavily managed catchments.

1 Introduction

Characterizing ecohydrological processes in sparsely monitored catchments with heterogeneous landscapes is inherently challenging due to spatially variable flow pathways and non-stationarity in rainfall inputs (Hrachowitz et al., 2013; McDonnell et al., 2007). This challenge can be even greater in catchments heavily modified by human activities, where a long and on-going history of disturbance can fundamentally alter processes and functioning (Marx et al., 2021). Distributed hydrological models are useful tools in addressing these challenges and are capable of capturing the dominant processes across spatio-temporal scales through regional parameterization (Fatichi et al., 2016). However, increasing model complexity to capture catchment heterogeneity makes it difficult to identify when models give “the right answer for the wrong reason” (Kirchner, 2006). In most catchments, rainfall and streamflow are the only available data for modelling. Streamflow-based calibration has therefore been the standard approach in hydrological modelling, leveraging the widespread availability of river discharge data to estimate model parameters across diverse catchments (Hrachowitz et al., 2013). However, calibration based on streamflow observations (single or multiple gauges) alone are usually insufficient to constrain hydrological model uncertainty, as certain parameters remain non-identifiable (Herrera et al., 2022). Consequently, simulations with multiple parameter sets can give equally plausible outputs, with process equifinality being a pervasive issue in model applications (Beven, 2006). Multi-criteria calibration, leveraging complementary datasets (e.g., soil moisture, ET, groundwater) in addition to streamflow, to mitigate this effect is increasingly common (He et al., 2019; Kuppel et al., 2018; Oliveira et al., 2021; Shah et al., 2021; Wu et al., 2023). However, classical hydrometric observations alone usually do not contain sufficient information to calibrate models and other observations such as isotopes are needed to improve understandings of storage and flux dynamics in catchments (Schilling et al., 2019).

Stable water isotopes (δ¹⁸O, δ²H) are proven tools that can help identify water sources, flow paths, and transit times, thus revealing process-based function in catchments (Holmes et al., 2023; Klaus and McDonnell, 2013; Sprenger et al., 2015). They are increasingly used as complementary datasets to streamflow, offering additional insights into catchment hydrological behavior that can aid model parametrization and calibration (Fenicia et al., 2008). In many cases, the integration of isotopes into the modelling process has advanced process representation, improving understandings of water partitioning and storage-flux interactions in heterogeneous landscapes (Birkel and Soulsby, 2015; Luo et al., 2024; McDonnell and Beven, 2014; Smith et al., 2022). Consequently, tracer-aided models (TAMs) have been increasingly applied worldwide (Jung et al., 2025).

Almost all TAM studies show inevitable trade-offs in model performance resulting from conflicting information in the streamflow and isotopes timeseries (Birkel and Soulsby, 2015; Scudeler et al., 2016; Wu et al., 2023). Such differences can highlight errors in model structure and inappropriate process conceptualization (Beven, 2006; McDonnell et al., 2007; Wu et al., 2025). In addition, using some observations as “soft data” (i.e. qualitative information or measured data that are not used in calibration) to constrain models can alleviate some of the above issues (Efstratiadis and Koutsoyiannis, 2010; Wu et al., 2023). However, failing to rigorously evaluate trade-offs between isotopes and streamflow in calibration risks producing structurally biased results, even if models achieve seemingly acceptable objective metrics for both datasets. Explorations of how these trade-offs in multi-criteria modelling using isotopes to help better understand hydrological processes and indicating further improve models are relatively rare.

Most TAMs focus on rural catchments (e.g. Soulsby et al., 2015) with limited anthropogenic disturbance (Yang et al., 2023), while complexities of ecohydrological processes are exacerbated in human-dominated systems where management measures can fundamentally alter hydrological connectivity and function (Wada et al., 2017). This creates critical unknowns in characterizing hydrological processes under anthropogenic impacts. In this regard, advancing tracer-aided methods to systematically evaluate hydrological dynamics at different spatial and temporal scales in heavily managed catchments can help improve modelling (Smith et al., 2021).

In this study, water stable isotopes (δ¹⁸O and δ²H) were used in a tracer-aided hydrological model (Spatially distributed Tracer-Aided Rainfall-Runoff, STARR), to help constrain estimates of ecohydrological partitioning and water balance compartments in four heavily modified sub-catchments of the Middle Spree catchment (MSC) in eastern Germany. These include the effects of agricultural irrigation, land use change, urbanisation and historic lignite mining with associated groundwater pumping to de-water coal seams. The area impacts a major national water resource, as the Spree river forms Berlin's water supply and ongoing pressures and intensifying climate change have the potential to threaten future water provision and ecosystem stability (Arndt and Heiland, 2024). Despite this significance, quantitative evaluation of ecohydrological processes in the MSC is currently limited, as records of intensive water use are not always available, historic impacts are often undocumented and parameterising these human influences in hydrological models are difficult. Therefore, this study aims to integrate streamflow and isotope calibration targets to provide a preliminary insight into ecohydrological couplings between storage and fluxes, as well as effects on the partitioning into runoff generation processes and ET fluxes in parts of the MSC. The study had the following specific objectives:

1. To assess the additional values of using seasonally-sampled isotopes in streamflow in addition to discharge in the simulation of heavily managed catchments.

2. To evaluate the contribution of isotopes in improving estimates of water flux partitioning and identifying the potential influence of unknown human management impacts.

3. Examining how management activities can bias ecohydrological models and advancing isotope-aided methods to disentangle process dynamics in human-dominated systems.

2 Materials and methods

2.1 Study catchment

The Mid-Spree catchment (MSC) is located in the SE (Southeast) of Brandenburg, Germany (Fig. 1). The 2806 km² sub-basin forms the middle part of the much larger Spree catchment (10 105 km²), accounting for 28.6 % of the entire catchment area. Within the MSC, the Spree River flows from Cottbus to Beeskow and through the Spreewald UNESCO Biosphere Reserve, which is an extensive wetland area. Climate is sub-continental with low precipitation and hot and dry summers (Pusch et al., 2009). Mean annual precipitation in the headwaters of the entire Spree catchment range from 600 to 1000 mm, decreasing to 556 mm in the MSC, making it one of the driest regions in Germany. Average monthly temperatures range between 19.3 °C in summer (June to August) and 1.9 °C in winter (December to February), respectively. Annual potential evapotranspiration, based on the FAO-56 Penman–Monteith equation (DWD, 2024), reaches 726 mm, making the MSC water-limited and highly susceptible to climate change.

Figure 1 (a) Elevation, river network and catchment borders, (b) land use, (c) geology, and (d) soil types of the four sub-catchments in the MSC.

The topography is flat and 80 % of MSC varies between 42 to 78 m a.s.l (meters above sea level), though the maximum elevation is 155.6 m a.s.l. The course of Spree River through the MSC has a very low gradient (0.027 %).

We selected four gauged sub-catchments in the MSC: the river Berste (gauged at Bruckendorf), the Wudritz catchment (gauged at Ragow), the Vetschauer Mill Creek (gauged at Vetschau), and the Dobra catchment (gauged at Boblitz), in the southern tributaries of the Spree River as the main study catchments (Fig. 1). They are mostly dominated by farming, particularly croplands (encompassing some of the most intensive farming areas in the MSC), with pasture and coniferous forests forming the other two major land uses in the four sub-catchments, accounting for 12 %–20 % and 18 %–30 % of the total area, respectively.

The geology of the sub-catchments is dominated by fluvial and meltwater sediments, especially along the river reaches, while artificial fill (from mining spoil) constitute a major part of the Wudlitz and Dobra catchments, as a result of historic lignite extraction (Landesregierung Brandenburg, 2024). Sandy soils (58.7 % of the area) and clay- sands (29.7 %) dominate these sub-catchments, while peat soils are sporadically distributed.

The four sub-catchments, like much of the MSC, were influenced by intense lignite mining activities between 1960–1990 (Pusch and Hoffmann, 2000). Pumped mine water (sump water), from de-watering former open-cast mines in the Altdoberm, Schlabendorf and Seese regions was discharged to the southern tributaries of the Spree River. After the sharp decline in lignite production in the early 1990s, discharge volumes from southern tributaries in the MSC to the Spreewald are believed to be close to the pre-mining situation since 2018 (Landesregierung Brandenburg, 2024). However, more generally in the Spree catchment, the lowered groundwater table leads now to a high risk of lower river flows and water shortages in the MSC, and constitutes a threat to Berlin's water supply further downstream (Arndt and Heiland, 2024). At present, drinking water in the studied region is extracted from groundwater wells. Agricultural and industrial water withdrawals are not permitted in the Vetschauer and Dobra catchments, while this extraction from the surface water for irrigation are limited to spring and summer at the Berste and upstream area of the Wudritz catchment. The Wudritz and Dobra catchments were most severely affected by historic mining activities several decades ago, and some parts of the mined areas have been converted into small lakes, with the Wudritz catchment having the largest surface water area among the four catchments (Table 1) (Landesregierung Brandenburg, 2024).

Table 1 Catchment characteristics of the MSC and its sub-catchments (Berste, Wudritz, Dobra, Vetschauer).

2.2 The STARR model and adaptations to simulate low-relief and different land use

The spatially distributed tracer-aided rainfall-runoff (STARR) model (van Huijgevoort et al., 2016) integrates the general conceptual structure of the HBV-light hydrological model (Lindström et al., 1997; Seibert and Vis, 2012) in a distributed, gridded framework that enables flux tracking of water and tracers through catchments. It is operated in the PCRaster Python framework (Karssenberg et al., 2010). The STARR model includes a module for tracking of stable water isotopes and tracer mixing (van Huijgevoort et al., 2016; Ala-Aho et al., 2017a). As a fully distributed model, hydrological fluxes are simulated in each grid cell based on a simple reservoir structure and water balance equations (Fig. S1 in the Supplement). The canopy, soil and groundwater storages are each conceptualised as a single layer, with throughfall, seepage and capillary flux serving as the flow paths connecting them. The model has been successfully applied in TAMs across a range of catchments at multiple scales (0.2–2500 km²) in contrasting environments (Ala-Aho et al., 2017b; Correa et al., 2020). The basic hydrological components and a brief summary of the STARR model are given in the Supplement (Sect. S1), while van Huijgevoort et al. (2016) and Dehaspe et al. (2018) provide more detailed descriptions.

Some of the assumptions of the STARR model, which was originally developed for mountainous catchments and subsequently adapted for tropical and cold regions, are not applicable in the lowland catchments studied here (the detailed equations used in modified STARR model are shown in Table S1 in the Supplement). For example, a topographical wetness index was not used to separate hillslope and lowland areas for runoff generation as in previous applications. Further, runoff routing was determined by the Manning equation (Chow, 1959), rather than assigning a pre-defined velocity.

Further, fractionation of stable water isotopes by evaporation in soil and interception storage was adapted to follow the Craig-Gordon model (Craig and Gordon, 1965), rather than being simulated by empirical representations (as in van Huijgevoort et al., 2016; Correa et al., 2020). The partitioning of evapotranspiration (ET) in the original STARR model was not applicable in the MSC as the isotopic composition of evaporation was sometimes more enriched than calculated isotopes of ambient atmospheric vapour (Chakraborty et al., 2018). Therefore, the partitioning method from HYDRUS-1D (Simůnek et al., 2013) was used here, where the transpiration originates from, and is linearly related to, soil and groundwater saturation. In order to keep parameter consistency, we adapted the interception module to the one developed in the HYDRUS-1D and the EcoPlot models (Stevenson et al., 2023), with the interception volume and the maximum capacity being controlled by LAI. Channel grid cells were distinct in representing different runoff and routing processes. Roughness (Manning coefficient) used in the kinematic wave equation was defined as 0.025 in the channel grid cells as recommended in Chow (1959), while other non-channel grid cells were assigned the values similar to van der Sande et al. (2003) and according to the land use (Table S2). The channel width used in the kinematic wave equation was estimated from recent Google earth maps (3 August 2024), while the width for non-channel grid cells was defined as the grid cell size (i.e., 500 m). Other parameters were the same for all grids. Explicit parameterization of anthropogenic factors (i.e., water withdrawls in the Berste and upper Wudritz, operation of restored lakes in Wudritz and Dobra) was excluded due to insufficient data, though the modelling results allowed these influences to be hypothesized.

2.3 Data acquisition

2.3.1 Forcing datasets

The spatial resolution of the model grid was defined at 500 m as a trade-off between adequate spatial detail and computation time (Smith et al., 2021). All datasets listed in Table 2 were downscaled or upscaled to the same resolution for consistency.

Table 2 Overview of the datasets used in this study.

The meteorological inputs were acquired from twenty weather stations in or near the study catchments – i.e., measuring precipitation (P), temperature (T), relative humidity (RH) – and grid products of potential evapotranspiration (PET) operated by the German Weather Service (DWD) were used, and the station datasets were linearly interpolated to spatially distributed (500 m) inputs. The 8 d composite LAI dataset at 500 m resolution  MCD15A3H V6.1) was used to characterise LAI dynamics. The dataset was accessed through Google Earth Engine (GEE). The cloud masking process was based on GEE and linear resampling at daily resolution was conducted, corresponding to the other input datasets. A global product to estimate the stable water isotopic composition (Interannual Monthly Mean values) of rainfall from (Bowen et al., 2005) was used and set as daily rainfall isotope input to the model constant in each month and equal to the monthly product value. This monthly value showed a consistent seasonal pattern with measurements from a rainfall station (daily resolution and aggregated to monthly values) downstream of the study area (Demnitzer Millcreek Catchment (DMC); 52°23^′ N, 14°15^′ E; Landgraf et al., 2022a), indicating the representativeness of the modelled rainfall isotope input (Fig. S2).

2.3.2 Datasets for model calibration and evaluation

Daily discharge from 2018 to 2023 and seasonal streamwater isotope data were collected during 2021–2023 in the outlets of the four catchments at the Bruckendorf (Berste), Ragow (Wudritz), Boblitz (Dobra), and Vetschau (Vetschauer mill creek) (Landesregierung Brandenburg, 2024). Gauging stations (Fig. 1) were used for model calibration (Table 2).

Stable isotopes were sampled (via instantaneously collected “grab” sampling) every season over 3 years (2021, 2022, 2023) at the river outlet of the four catchments, and we used cavity ringdown spectroscopy (Picarro L-2130i, Picarro Inc., CA, USA) to analyse the ratio of the heavy stable isotopes of oxygen (¹⁸O) and hydrogen (²H) relative to the Vienna Standard Mean Ocean (VSMOW) (for detailed sampling and pre-process procedure refer to Chen et al., 2023).

Further, the MODIS ET (MOD16A2GF from GEE) and PML ET (Zhang et al., 2019) 8 d composite products were compared with simulation results to evaluate evapotranspiration simulations, and qualitatively check the model performance. Despite uncertainties, the PML product aligned better with flux tower (51.8922° N, 14.0337° E) records in the Spreewald (Table 2) indicating its usefulness as a comparator for modelling results (Fig. S3). Groundwater (GW) levels from 3 wells distributed in the studied region were collected (Landesregierung Brandenburg, 2024) and the annually varying pattern of each well were used to further identify the robustness of the model.

2.4 Model parameterisation, sensitivity analysis and calibration

2.4.1 Model parameterisation

The number of calibrated model parameters was minimised and therefore some of the parameters were assigned fixed values (Table S2). In total, 35 parameters were included for calibration and the assigned ranges for calibration were mostly adapted from previous applications of the STARR model, with some adjustments appropriate to the characteristics of the study catchments (Table S2). Parameters representing soil characteristics were distributed according to land use types (i.e., urban, water, pasture, cropland, shrub, forest) given the close correlation between land cover and soil type in the region (see Smith et al., 2021). The fluxes from interception to soil storage (throughfall and stemflow) are different under forest and non-forest land use, as contrasting canopy characteristics affect the rainfall partitioning (Guevara-Escobar et al., 2007), and the corresponding parameters were determined separately.

2.4.2 Sensitivity analysis (SA)

We employed the Morris Method (Morris, 1991) to identify the most sensitive parameters. The calculation used the SAFE tool (Sensitivity Analysis for Everybody, Pianosi et al., 2015). The elementary effects (reflecting the sensitivity of each parameter) were calculated through perturbing the starting parameter by a certain variation based on a radial one-at-a-time strategy. Kling–Gupta efficiency (KGE) (Gupta et al., 2009) of simulated streamflow or stream water isotopes, corresponding to calibration scheme, was used as the objective function in calculating the elementary effects (Table 3). The Latin-Hypercube sampling method (Pianosi et al., 2015) was selected to determine the starting parameter and following perturbation. The mean elementary effect of each parameter was used to indicate the sensitivity of the corresponding parameters. Parameters related to pastures, croplands, and forests were more sensitive than others, as they cover majority of the catchment area (Table 4), and the estimated ecohydrological fluxes are mainly shown for these three land uses.

Table 3 Five calibration schemes based on measured streamflow and isotope at the outlet of the four catchments. Sensitivity analysis was conducted separately on isotope and discharge measures in schemes 2–5.

Table 4 Performance metric values (KGE) of simulated spatially averaged ET in the scheme 1 and balanced solutions of schemes 2–5 based on remote sensing products at the studied catchments.

2.4.3 Model calibration

The modified STARR model was run at daily time steps for the period from 2014 to 2023. A 4-year spin-up period (2014–2017) was applied, and the remaining 6 years were used for calibration. The multi-objective non-dominated sorting genetic optimization algorithm II (NSGA-II) (Blank and Deb, 2020; Deb et al., 2002) was applied to derive the Pareto-optimal solutions. Five distinct calibration schemes were conducted based on measurements (i.e., streamflow and/or stream isotopes) at the outlet of the four catchments and are detailed in Table 3. The first scheme (discharge-only based calibration) was calibrated only on the streamflow of all four catchments (multi-gauged streamflow calibration) using the KGE values of simulated streamflow as the objective function. Considering the potential for heterogeneity in the hydrological functioning of each catchment and to assess the additional information content of the isotope data in the modelling, calibrations (the other four schemes as in Table 3) were also carried out in each catchment independently (Scheme 2: Berste; Scheme 3: Wudritz; Scheme 4: Vetschauer; Scheme 5: Dobra) and based on KGE values of simulated streamflow and isotope (isotope-aided cal generations and 280 individuals in the population per generation, and the first pareto front in the 80th generation were employed as the solution. Five independent calibrations were conducted for each calibration scheme to minimize potential biases. Balanced, discharge-dominated, and isotope-dominated solutions were defined in the isotope-aided calibrations (schemes 2–5), five parameter sets from each of the five independent calibrations (25 parameter sets in total) were selected for each of these three solution types. These three solution types corresponded to parameter sets in the middle – with minimum Euclidean distance to the best KGE of discharge and isotopes, i.e., (KGE_discharge = 1.0, KGE_isotope = 1.0) – and the edges of the Pareto front, respectively. In the scheme 1, only balanced solutions (minimum Euclidean distance to the best KGE point, i.e., KGE_Berste = 1.0, KGE_Wudriz = 1.0, KGE_Vetschauer = 1.0, KGE_Dobra = 1.0) were employed. The aggregated 25 parameter sets for each solution type (i.e., balanced, discharge-dominated, isotope-dominated) encompass a broad range of parameter variability and allow a robust comparison among solution types. No validation period was employed, as the present study aimed to better constrain ecohydrological processes in the MSC, rather than for forecasting applications. Also, others have argued that calibrating to the full available data and skipping model validation can be more robust than the traditional split-sample approach in hydrological modelling (Shen et al., 2022).

3 Results

3.1 Model performance

3.1.1 Discharge simulations

Discharge-only based calibrations (scheme 1) and the balanced solutions of the isotope-aided calibrations (schemes 2–5) successfully captured discharge dynamics in each catchment with average KGE > 0.5 (Fig. 2). This performance was better (KGE > ∼0.8) in discharge-dominated solutions of the isotope-aided calibrations (Fig. S4). Contrasting flow regimes were shown in each catchment, e.g., Berste presented low base flows and Vetschauer had quicker responses to rainfall. Therefore, scheme 1 underestimated peak flow, showing the trade-offs in balancing performance across the four catchments. Trade-offs between discharge and isotopes objectives within each catchment (schemes 2–5) also led overestimated base flow (particularly in Berste and Wudritz) in the balanced solutions of schemes 2–5 (Fig. 2).

Figure 2 Discharge simulations at the outlet of each catchment based on the balanced solutions in the calibration scheme 1 (first panel a1, b1, c1, d1) and schemes 2–5 (second panel a2, b2, c2, d2). (a) Berste; (b) Wudritz; (c) Vetschauer; (d) Dobra.

3.1.2 Isotope dynamics

The stream isotope signatures in the four catchments showed contrasting characteristics (Fig. 3). Overall, apart from the Berste, streamwater isotopes in each catchment plotted below the local meteoric water line (LMWL), reflecting differences in fractionation. The similar alignment of isotopes along a shared local evaporation line indicates comparable atmospheric moisture demand among the catchments (Fig. 3). The Berste exhibited the most depleted isotopic signature (δ²H: $\sim -\mathrm{62}$ ‰), while Wudritz was the most enriched catchment (δ²H: $\sim -\mathrm{38}$ ‰), and their enrichments positively correlated with the extent of surface water area, although all catchments have low surface water coverage (<8 %) (Fig. 3). Simulations in all solutions of each scheme reproduced the seasonal isotope dynamics, with summer enrichment and winter depletion (Figs. 3 and S5). However, the variability of isotopes was overestimated in scheme 1 and in the discharge-dominated solutions of schemes 2–5, although mean simulations and measurements were comparable and much better performance were shown in the latter. This showed that different flow paths, mixing processes and fractionation effects in the catchments were likely unrealistic in the discharge-only based calibrations (scheme 1) or discharge-dominated solutions in isotope-aided calibrations (schemes 2–5). The balanced solutions or isotope-dominated solutions in schemes 2–5 yielded much more consistent simulations of the isotope dynamics, although this came at the cost of much poorer discharge performance (Figs. 3 and S5).

Figure 3 Stream water isotope simulations at the outlet of each catchment based on the balanced solutions in the calibration scheme 1 (first panel a1, b1, c1, d1) and schemes  –5 (second panel a2, b2, c2, d2). (a) Berste; (b) Wudritz; (c) Vetschauer; (d) Dobra; and (e) Dual isotope from January 2021 to December 2023 in each catchment.

3.2 Trade-offs between streamflow and isotope-based calibrations

3.2.1 Pareto front of simulations

The Pareto fronts of simulations from schemes 2–5 show the range of potential solutions of streamflow and isotopes (Fig. 4). Berste resulted in lower KGE of both isotopes and streamflow than the other three catchments, while the balanced trade-off solutions in Vetschauer were the most robust (Fig. 4). On each Pareto front, increased KGE of simulated discharge was accompanied by a slight decrease in the KGE of simulated isotopes, until the discharge KGE reached a certain threshold, beyond which the isotope KGE dropped sharply. This threshold was approximately 0.7–0.8 at Wudritz, Vetschauer and Dobra, but only around 0.5 at Berste. This unsatisfactory performance at Berste reflected the stronger conflicting information provided by isotope and streamflow data under the existing model structure and potentially reflected unknown impacts resulting from management (see Discussion below).

Figure 4 Pareto fronts of simulations from all five independent calibrations in schemes 2–5.

3.2.2 Impact of calibration on contributing runoff generation processes

Runoff generation from soil storage was the most sensitive hydrological processes in the simulations in all calibration schemes, while groundwater (GW) contributions were highlighted by isotopes in Berste and Vetschauer (Fig. 5). Contrasting GW contributions to discharge were shown in different solution types of the isotope-aided calibrations (Fig. 6). In the discharge-dominated solutions of schemes 2–5, runoff was mainly sourced from the shallower soil stores (soil storage in the STARR), while this contribution decreased in the isotope-dominated solutions with groundwater sources being more important (Fig. 6). Summer peaks of groundwater contribution shown in the isotope-dominated and balanced solutions aligned better with the known hydrological regime of the local region (e.g. Ying et al., 2024), compared to discharge-dominated solutions (Fig. 6). Among the four catchments, higher groundwater contributions were shown in Berste and Wudritz, reaching 60 % in isotope-dominated or balanced solutions (Fig. 6).

Figure 5 Sensitive parameters in each calibration scheme (the top six sensitive parameters are shown shaded). I and Q in the y axis represent isotopes and discharge; (p), (c), (f) represent pasture, croplands, forest, respectively. kS, kG LP, α, K control processes of runoff from soil, runoff from groundwater, Actual ET, interception and ET partitioning, respectively; and FC is the soil water storage capacity.

Figure 6 The mean ratio of seasonal runoff from groundwater storage at each catchment (winter, spring, summer, fall of each year along the x axis). The results were calculated from the discharge dominated, balanced and isotope dominated solutions in schemes 2–5.

3.2.3 Impact of calibration on ET estimates and ET partitioning

The parameters controlling the return flux of water to the atmosphere (ET) were sensitive to the use of both discharge and isotopes in calibration (Fig. 5). Further, spatially averaged ET in the balanced solutions of all schemes aligned well with both remote sensing products (MODIS and PML), supporting the robustness of the calibrations, although the trade-offs between discharge and isotopes at Berste resulted in a degraded performance in the Scheme 2 (Table 4). All solution types in the isotope-aided calibrations generally showed only slight differences, except for croplands in Berste and Wudritz, where the discharge-dominated solutions simulated higher ET (Fig. S6). ET showed strong seasonality with peaks in spring and summer in all solution types (Figs. 7 and S6). This was consistent with the satellite ET observations. However, balanced solutions showed higher springtime ET yet lower modelled ET in summer across all catchments compared to the remote sensed products, due to limited soil moisture. These differences decreased during wet summers (e.g., 2021 and 2023) (Figs. 6 and S4).

Figure 7 Seasonal total ET and transpiration ratio (T $/$ ET) from the balanced solutions at (a) Berste in scheme 2; (b) Wudritz in scheme 3; (c) Vetschauer in scheme 4; (d) Dobra in scheme 5, and RS products. Bars and scatters indicate ET and T ratios, respectively. Each season contains three bars, ordered from left to right as follows: balanced solutions, MODIS and PML ET products.

ET in the discharge-dominated solutions of schemes 2–5 presented a relatively uniform spatial distribution, while the balanced or isotope-dominated solutions retained more pronounced spatial variability, aligned with the MODIS estimates (Figs. 8 and S7). However, no solutions successfully captured the distributions of ET in the MODIS or PML. In addition, high ET in the PML product was more consistent with croplands distribution, which indicated the potential influence of irrigation in croplands and were not shown in the MODIS, reflecting the uncertainty of accuracy of spatial distributions in these remote sensing products (Fig. 8).

Figure 8 Spatial distribution of ET in the four catchments. (a) Balanced solutions in schemes 2–5; (b) MODIS ET; (c) PML ET. Suffix “1” and “2” in titles of each subplot mean dry year of 2019 and wet year of 2023, respectively.

ET partitioning (into E and T) became more sensitive for simulations constrained by isotopes as the corresponding parameters (K) ranked higher compared to Scheme 1 or discharge-dominated solutions in Schemes 2–5. This is consistent with the importance of fractionation in the variation of isotopes (Fig. 5). In all catchments, the transpiration ratio from each land use had similar temporal dynamics (summer peaks and winter troughs), and was higher in forest than in the other two land uses (Figs. 7 and S8). Berste, Dobra and Vetschauer yielded higher transpiration ratios, with 70 %–80 % and 40 %–60 % in summer and winter, respectively, compared to Wudritz. Further, Wudritz was ∼10 % lower (compared to other catchments) in transpiration ratio, aligned well with the PML estimates. In addition, Berste and Wudritz showed contrasting performances among different type of solutions (discharge-dominated, balanced or isotope-dominated solutions), with discharge-dominated solutions showing higher ratios (Fig. S8).

3.3 Quantification of water balance components

The partitioning into different water balance components was relatively consistent across all catchments, irrespective of the calibration schemes (Fig. 9). ET was the dominant output flux, especially during the spring and summer period (>90 %). Simulations showed that soil storage played the major role in supplying water for ET in spring and early summer and was subsequently replenished during autumn and winter. Discharge accounted for a minor proportion of the water balance (∼5 %). Variations of groundwater storage in each season were small (<1 % of the total water balance) in simulations calibrated on discharge alone, while it increased to ∼5 %–30 % (catchment dependent) when isotopes were included in calibrations (Fig. 9). A clear seasonal pattern of groundwater dynamics was shown in Vetschauer, with declining storage during the spring and summer and replenishment in winter and fall. In contrast, the groundwater storage kept increasing in each season at Berste, while this storage only occasionally decreased in summer at Wudritz and Dobra (Fig. 9). The annual increased in subsurface storage aligned with the gradually recovery of the groundwater level in this region following the 2018–2019 drought (Fig. S9).

Figure 9 Seasonal water balances at (a) Berste; (b) Wudritz; (c) Vetschauer; (d) Dobra, based on the balanced solutions in the calibration scheme 1 (first panel a1, b1, c1, d1) and schemes 2–5 (second panel a2, b2, c2, d2). Positive bars represent water sources, while negative bars are water losses. “Channel storage” and “Interception storage” were too small compared to other components to show in the figure.

4 Discussion

The key novelty of this study lies in evaluating the use of multi-year, seasonally sampled isotopes for calibrating a TAM to ET-dominated and heavily managed catchments. Previous catchment-scale TAM applications have generally been in smaller experimental catchments with limited human influences (Nan and Tian, 2024). Whilst such work has demonstrated the value of isotopes in constraining hydrological models, the present study has shown similar potential benefits in heavily managed regions. Furthermore, the trade-offs between using discharge and/or isotopes in model calibration showed that isotopic information can still provide valuable preliminary process insights. The observed inter-catchment differences in isotope variations were influenced by the interaction of dominant runoff sources, associated mixing processes with water storages and the effects of evaporative fractionation; capturing these in simulations were, in turn, linked to modelled discharge performance. The response of runoff to rainfall, seasonal flow variations and baseflow contribution derived from discharge-calibrated models showed inconsistency with isotope-inferred processes. Given that contrasting trade-offs were identified among catchments within the same region, these discrepancies are likely to in part reflect undocumented influences of human management and provide a potential pathway to model improvement.

4.1 Benefits of incorporating isotopes into model calibration

Streamflow is most widely used in hydrological model calibration due to its ready availability. Here, similar to many previous studies, streamflow dynamics and overall runoff volumes were effectively captured by discharge-only based calibrations despite the ET-dominance. However, slight overestimations of flow occurred during winter peaks, with underestimations during summer low-flows. In addition, simulated spatially averaged ET (weekly flux) aligned closely with RS ET products (e.g., MODIS or PML), with KGE > 0.5. Although groundwater recharge was implausibly low, approximate partitioning of precipitation into green and blue water fluxes was accurately constrained by streamflow calibrations. Similarly, in some other literature, multi-criteria calibrations incorporating both streamflow and RS ET demonstrated marginal improvements in annual ET or streamflow simulations, compared to discharge-only based calibrations, albeit with reduced parameter equifinality (Oliveira et al., 2021; Shah et al., 2021), although in other cases they have not been beneficial (Ala-Aho et al., 2017b).

However, previous work using tracers has shown that based on calibration on streamflow alone can result in the wrong conceptualization of hydrological processes (Ala-Aho et al., 2017a; van Huijgevoort et al., 2016; McDonnell and Beven 2014). Discharge-only based calibrations normally exhibit significant uncertainty in representing internal hydrological processes, owing to high degrees of model freedom and non-identifiable parameters (Herrera et al., 2022). This was evidenced by poorly simulated isotopes under scheme 1. Moreover, catchment discharge, as an output integrating upstream hydrological processes, provides only limited insight into spatially distributed partitioning of water volumes and flow paths (Fatichi et al., 2016). This limitation becomes more pronounced as models incorporate more complex spatial disaggregation and physics-based process representations (Sun et al., 2017). A multi-gauge regional streamflow calibration, as used in scheme 1, has been shown to improve the spatial representativeness of simulations and provides a greater constraint on relationships between catchment characteristics and streamflow dynamics (Liu et al., 2020). However, internal mixing and runoff generation processes (evidenced by isotope dynamics) under scheme 1 did not show much difference to the single-gauge calibrations (i.e., discharge-dominated solutions in schemes 2–5). Improvements on simulations of subsurface processes were also not clear in other literature (Wanders et al., 2014). Furthermore, sensitivity analysis revealed that discharge was more sensitive to soil drainage process in all three major land uses (i.e., forest, croplands, pasture) of the four studied catchments, and parameters controlling other hydrological processes were less identifiable.

Incorporating stable water isotopes improved modelling of ecohydrological partitioning and storage dynamics in our study. Multi-criteria calibration using both streamflow and isotopes enables a wider exploration of the parameter space to satisfy multiple calibration objectives (Holmes et al., 2020). Subsurface processes in the deeper layer (e.g., discharge from groundwater in the present modelling) are inherently challenging to constrain using near-surface observations (e.g., ET, soil moisture), and these processes are usually poorly understood (Beven, 2006). Stable water isotopes are powerful in this regard as they integrate the cumulative effects of water flow paths and mixing across all hydrological storages (Godsey et al., 2010; Tafvizi et al., 2024). Our TAM revealed that discharge-only based calibration or discharge-dominated solutions in isotope-aided calibrations failed to capture such mixing processes, as their simulated isotopes deviated far from observations. Modelling results that reconciled simulated and observed isotopes likely better represented mixing between soil water and groundwater storage. This generally increased runoff generation from groundwater, characterized by more depleted isotopes, and modulated simulated outflow signatures, flattening the seasonal isotopic variability, and thus, giving more plausible process representation. This aligns with previous findings where dominant subsurface flows produce subtle isotopic variations (Iorgulescu et al., 2007; Oerter and Bowen, 2019) and where contrasting isotope signatures across water stores are key to disentangling water mixing processes (Kirchner, 2003).

The fractionation of stable water isotopes is governed by evaporation, and measured isotopes at the outlet reflected the aggregated influence of evaporation across the whole catchment. Parameters controlling ET processes (e.g., FC and LP) exhibited similar sensitivity to isotopes and discharge, supporting the value of isotopes as a constraint on ET dynamics in simulations. This is similar to the incorporation of remote sensed ET in calibrations found in other literatures (Oliveira et al., 2021). In addition, total seasonal ET volumes across the catchments showed minimal variations among different solution types in the Pareto front, implying that streamflow and isotopes similarly constrain bulk ET estimates. The use of relatively coarse isotopic datasets and inter-annually averaged rainfall isotope inputs in this study likely introduced some uncertainty in ET process constraints, and higher-resolution data would better capture temporal variability in fractionation intensity (Sprenger et al., 2017). Nevertheless, the results are reasonable and consistent with broader findings in the literature: while isotopic tracers are widely recognized for improving estimates of $E/T$ partitioning (Gibson and Edwards, 2002), their utility in refining total ET quantification prior to discharge remains less clearly demonstrated, particularly at larger spatial scales.

Despite this, isotopes are still likely to provide more additional information than ET observations for constraining ET dynamics at the catchment scale. Although calibrations based on discharge and temporal patterns of spatially aggregated ET can help constrain total ET in an efficient manner (Oliveira et al., 2021; Shah et al., 2021); and spatial ET patterns have also helped improve spatial representation of ET (Dembélé et al., 2020), resolving the underlying components contributing to ET with contrasting residence times (e.g. interception, soil evaporation, transpiration etc.) remains difficult. In this regard, isotopes offer additional constraints separating the sources of individual hydrological components, as well as the capability of tracking sources of ET, as has been demonstrated in numerous previous studies (Martinelli et al., 1996; Ren et al., 2024). In the present study, the dominant contribution of groundwater to streamflow inferred from the isotope-aided calibrations indicated reduced storage in the shallow subsurface layer (soil storage in STARR), which in turn leads to improved estimates of total ET and ET water sources compared to discharge-only calibrations.

4.2 Isotopic insights into water flux partitioning in managed catchments

In the studied catchments, runoff was predominantly generated from groundwater and soil storage, reflecting a subsurface-process dominated flow regime. This pattern aligns with observations across much of the Spree catchment, where subsurface-driven runoff mechanisms are widespread (Chen et al., 2023). Runoff in the upper Spree catchment is predominantly groundwater-driven (Kröcher et al., 2025), a pattern consistent with the simulated performance in the studied region after incorporation of isotopes. In addition, the influence of groundwater in the studied four catchments may potentially increase, as the opencast mines have been closed for 30 years and it is unclear to whether groundwater levels have fully stabilized and how they affect each of the study catchments (Kröcher et al., 2025). The increased temperature and shift of precipitation from summer to winter due to climate change also possibly leads to increased ET during winter and spring reducing discharge and groundwater recharge accordingly. This could result in contrasting water distributions in each season, and intensifying negative climatic water balance in the local environment (Pohle et al., 2012).

The transpiration ratios in the studied region were roughly 65 % in the PML estimates, which is similar to other studies in the Spree catchment (Landgraf et al., 2022b). However, simulated transpiration ratios at Berste, Vetschauer and Dobra were ∼10 % higher than the PML estimates, consistent with the limited fractionation shown in the streamwater isotopes. The temporally coarse sampling potentially underestimated the seasonally aggregated fractionation, but the general seasonal variations were captured in the data. In addition, uncertainties in this remote sensing product cannot be discounted. In contrast at Wudritz, the simulated ratios were more consistent with the PML values. The mismatch in transpiration ratios may have resulted by underestimated evaporation from calibration with seasonally sampled isotopes, though this may be counterbalanced by unaccounted surface water evaporation from lakes in the former mining area. While evaporation from small open water bodies has negligible impacts on overall the catchment water balance, it can play an important role in isotopic enrichment (Birkel et al., 2011). For instance: in the Wudritz (7.7 % surface water area), isotopic concentrations in simulations were likely underestimated due to unmodeled surface water evaporation. The non-linear relationship between evaporation rate and isotopic enrichment, as described by the Craig–Gordon model (Craig and Gordon, 1965), explains this dynamic: early-stage evaporation induces stronger isotopic enrichment, approaching a threshold under constant environmental conditions (e.g., humidity, temperature). Thus, even relatively minor surface water evaporation can bias isotopic signatures which then impacts ET partitioning simulations.

Further, the small weight of isotopes in a scalar function combining multiple objectives meant that it was possible to help disentangle ET processes (Wu et al., 2023). This was also illustrated by the minor adaptations along the Pareto front at Vetschauer and Dobra. The normally contrasting transpiration ratios between different land uses (Schlesinger and Jasechko, 2014) were consistent with our simulations, and divergent LAIs in each land use explains these differences (Cao et al., 2022), although the PML estimation presented similar transpiration ratios among each land use. However, the transpiration ratios could be underestimated in specific land uses due to unparameterized irrigation effects (Paul-Limoges et al., 2022), and the local water use efficiency needs to be further evaluated to improve future modelling.

Despite uncertainties introduced by relatively coarse sampling, isotopes were still valuable in water partitioning (i.e., ET and runoff) in such heavily managed catchments, which has been shown to be an advantage compared to other potential calibration targets (e.g., ET, soil moisture; Wu et al., 2023). The preliminary assessment of ET partitioning through isotope-aided modelling provides evidence for assessing the water-use efficiency of different land uses, especially croplands. Isotope data at finer temporal resolution would likely help better constrain ET partitioning in such heavily modified catchments, as they do in less disturbed environments (Tunaley et al., 2017; Soulsby et al., 2015). In addition, precipitation isotopes were from a global data product (Bowen et al., 2005). The use of interannually averaged monthly values likely failed to capture short-term climate variability or anthropogenic influences (e.g., fossil fuel-derived vapor), and may introduce uncertainties in water source apportionment (Xia et al., 2024; Yang and Yoshimura, 2024), though given the groundwater dominance and older (>5 years) ages of runoff (Smith et al., 2021), the effects are likely to be small. Again, high resolution local data would be advantageous to improve such estimates.

4.3 Insights from contrasting discharge-isotope trade-offs among the four catchments

Trade-offs in the calibrations between streamflow and isotopes are a well-known feature of TAM (Holmes et al., 2020), though their severity varies across applications. The extent of these compromises depends on the model's structural flexibility to assimilate additional constraints (Holmes et al., 2020) and the weightings given (Wu et al., 2025). Whilst some applications based on spatially-distributed models presented relatively minor conflicts in the information content of different calibration targets (Kuppel et al., 2018), significant degradation in streamflow performance has been found in lumped models after incorporating isotopes in calibration (Fenicia et al., 2008). In our study, the lack of explicit representation of some anthropogenic drivers (e.g., water abstraction, irrigation etc.) in the model emerged as a likely contributor to trade-offs. Without explicit parameterization for these factors, the model often failed to fully capture flow dynamics in the catchments.

There was evidence for such systematic biases. For example, discharge-dominated solutions in isotope-aided calibrations produced a soil water-driven runoff regime that generally underestimated baseflow. The STARR model derived a fast draining process in the shallow subsurface layer (soil storage in the STARR) to simulate flattened isotope variations. Water in the shallow subsurface flows to the outlet through direct discharge (faster draining lead to lower evaporation and lower fractionation) and to the groundwater through seepage (high groundwater storage flattening isotope variations). In terms of discharge, high groundwater storage would typically result in flattened seasonal variations and elevated baseflows, whereas rapid draining of soil storage to discharge resulted in quick runoff response to rainfall, as illustrated by the isotope-dominated solutions (Fig. S4). However, in the studied region, only the Vetschauer catchment showed both a rapid response to rainfall and high baseflow, and therefore showed relatively minor trade-offs between isotopes and discharge in the modelling. In contrast at the Berste catchment, discharge exhibited low base flow and pronounced seasonal variations but lacked a rapid response to rainfall events, which was inconsistent with the hydrological processes inferred by isotope observations. As all four catchments are located within the same region, with broadly similar soil cover, it might be expected that the soil hydrological response would be similar among the four catchments. However, the Berste, with the largest area, had the lowest base flows and slower responses to rainfall events compared. Considering the Berste catchment has the most extensive cropland coverage and high (but unknown) irrigation withdrawals, these trade-offs could be potentially explained by management measures. Moreover, the higher groundwater contribution at Berste and Wudritz, compared to Vetschauer and Dobra, might be explained irrigation water withdrawls. Extraction of groundwater to the surface could dampen isotope variations in stream water, although much of the water volumes are expected to evaporate during summer period. Therefore, the un-represented water extraction process likely introduces more conflicts between discharge and isotope observations.

These trade-offs would likely be mitigated by more process-based conceptualizations of anthropogenic impacts (e.g. irrigation), though at the expense of greater parameterization and higher uncertainty. The STARR model's simplified routing structure, where runoff from contrasting sources follows identical pathways to the outlet, diverges from more physics-based frameworks that better spatially differentiate the timing of contributions of overland flow, unsaturated zone flow, and groundwater flow. This conceptual routing can capture flow dynamics in the stream, but likely mis-represent flow partitioning and creating conflicts in multi-criteria calibrations (McDonnell and Beven 2014). An additional point is that STARR did not include any consideration of channel leakage and parameterization of the hyporheic flow could potentially and partly explain the flattened isotope variation and over-estimated baseflow simulations. This is a particular problem in former mining areas where recovery after historic groundwater is incomplete or has altered local groundwater surface water interactions increasing the possibility of streams having “losing” reaches with vertical leakage into groundwater (Pohle et al., 2025).

Further, STARR calculates ET solely based on water availability and partitions PET into PE and PT using a simplified single calibrated parameter. In contrast, more physically-based models, such as EcH2O-iso (Kuppel et al., 2018), explicitly simulate runoff generation processes and compute ET through coupled energy and water balance equations, and thus likely better integrate isotopic information. While such enhancements improve a model's ability to simultaneously reconcile flow velocity and celerity, they also introduce greater parameter uncertainty (Herrera et al., 2022), and a trade-off between benefits of increased model complexity and increased simulation uncertainty will result (Schilling et al., 2019). Simple explicit modeling of water extraction for irrigation (e.g., channel-to-cropland transfers) could alleviate trade-offs observed in this study. However, sparse records of withdrawal volumes, associated irrigation patterns and potential unauthorized practices, as well as political sensitivities in data sharing may limit practical modelling.

Despite having relatively limited streamwater isotope samples compared to some other studies, their availability was invaluable for this intercomparison and revealed some significant difference between catchment function (Fig. 3). Since the evaluation of trade-offs between discharge and isotopes primarily relied on the seasonal patterns of isotopes, the 3 years of seasonally sampled isotopes used for calibration and the monthly rainfall isotope inputs provide a sufficient basis for evaluation for the trade-offs. In addition, the potential anthropogenic impacts noted above were identified through the comparisons among the four catchments in the same region, further supporting the value of the analysis. This evaluation of trade-offs between discharge and isotopes helped identify potential processes not adequately captured through discharge alone in such a heavily managed catchment.

5 Conclusions

Stable water isotopes are valuable tracers for identifying water sources, offering the potential to constrain ecohydrological models. While model calibration in natural catchments typically exhibits slight trade-offs between isotopic signatures and conventional hydrological variables (usually discharge), this study advances a novel perspective on the benefits and challenges of integrating isotopes in catchments heavily impacted by human activities. Using the conceptual-based, fully-distributed TAM in the STARR model, we calibrated both isotopes and streamflow without explicitly parameterizing anthropogenic disturbances to investigate three critical issues: (1) the potential of using discharge alone in calibration to mislead process interrelations from simulations under anthropogenic stress, and (2) the adaptability and value of seasonally sampled isotopes in such contexts, and (3) the likely impacts of human interventions on model performance. We studied four sub-catchments of the Middle Spree (Berste, Wudritz, Vetschauer, and Dobra), subjected to contrasting anthropogenic pressures (long-term mining impacts, seasonal water withdrawals), and derived Pareto-optimal solutions to help disentangle the additional insights provided by isotope-aided calibration compared with streamflow alone.

The results demonstrate that discharge-only based calibrations could mis-represent runoff partitioning processes, especially in catchments with water withdrawals for irrigation, while isotopes helped identify implausible simulations by being more informative on process representation. For example, higher groundwater contributions were identified through isotope-aided calibrations. Isotopes also helped to disentangle ET partitioning, although seasonally sampled isotopes could underestimate seasonally aggregated fractionations, and the absence of explicitly parameterizing the hydrology of small mining lakes could potentially mislead the results. Strong trade-offs between isotopes and streamflow in calibrations arise in such anthropogenically-impacted catchments, where unquantified epistemic errors in streamflow observations caused by human activities can compromise model reliability.

This study highlights how unaccounted anthropogenic activities can alter model interpretations and underscores the complementary role of isotopes and TAMs in refining simulations under complex human-environment interactions, although only seasonally sampled isotopes were employed. While distinct trade-offs between streamflow and isotopes were observed in the study catchments, these simulations still offer informative insights into ecohydrological dynamics and partitioning in heavily impacted catchments, even when quantitative process evaluation remains challenging. In catchments subject to intensive anthropogenic interventions (e.g., altered water distribution via irrigation or withdrawals), the severity of streamflow-isotope conflicts and compromises in TAM may serve as an indirect diagnostic of human impacts on water partitioning. Representing anthropogenic effects in ecohydrological models is inherently difficult, particularly when historical data on water use or management practices are sparse. However, we demonstrated here that TAMs have valuable potential in such applications and encourage greater isotope sampling in anthropogenically influenced catchments.