The use of a shuffled complex evolution algorithm for calibration of a dynamic leaf emergence model leads to better spatial ETa predictions

Tree phenology as the study of the timing of seasonal biological events in trees, e.g. leaf development, flowering or leaf fall, has a crucial role in determining ecosystem processes and functions. For example, leaf emergence leads to increased evapotranspiration rates and thus water extraction from the soil by the plant (Leuning et al. 1991). The development of fruits is relevant for the reproduction of plants, but also for anthropogenic food production on agricultural land. Leaf fall, in turn, enables the remineralisation of biomass, but also acts as protection against frost (Esteso-Martínez & Gil-Pelegrín 2004). The influence of phenology is therefore not limited to the plant itself but extends far beyond that into diverse areas of the ecosystem. The phenology is influenced by a wide variety of environmental factors. For example, the sprouting of plants is highly sensitive to the amount of heat available, i.e. the temperature development in spring (Milthorpe 1959; Polgar & Primack 2011). The photoperiod also seems to influence phenological processes such as bud bursting in the genus Fagus L. or Aesculus L. and leaf primordia after late dormancy (Zohner & Renner 2015). In addition to the temperature, the precipitation and thus the moisture content of the soil are also an influencing factor for sprouting, for example in grassland (Shen et al. 2011). Studies have shown that phenological processes are good indicators of climatic variability, including those induced by climate change. This is due to the strong sensitivity of plants to environmental drivers, as evidenced by increasingly earlier leaf emergence and average flowering date (Bradley et al. 1999; Lu et al. 2006; Beaubien & Hamann 2011). Other indicators are their possible delay due to increased difficulties in reaching required chilling units, e.g. in Malus domesticus Borkh. or other fruit trees (Funes et al. 2016).

An investigation by Zhao et al. (2013) into the existing approaches for the simulation of dynamic phenology showed that there are mainly three different archetypes of plant phenology models: statistical, mechanistic and theoretical. Although they are all based on common driving parameters such as temperature, light availability or radiation, the model assumptions regarding the reaction behaviour of the plant are different. The statistical models mainly use heat units, photo-heat units or hydro-heat units as defining parameters. The regression models are often used to determine threshold values. Mechanistic models are usually based on forcing and chilling units, as well as corresponding modifications of these, as determining parameters for sprouting (Kramer 1994). Models based on these are, for example, the deepening rest model (Hänninen 1990) or the North Carolina Model (Clements & Ditommaso 2011). There are also theoretical models that are based less on forcing parameters and instead use aspects such as cost–benefit trade-offs or evolutionary relationships as a basis for explanation (Zhao et al. 2013). Models of this type are, for example, hormone and interaction-based models like the promoter–inhibitor model (Schaber & Badeck 2003) or genetic behaviour-based models like the phenology model developed by Wilczek et al. (2009). However, the majority of these models are mathematical–conceptual models that stand alone and are often not integrated into complex water balance models. This leads to only a few existing water balance models that can simulate dynamic phenology (Kochendorfer & Ramírez 2010; Dolschak et al. 2019; Chen et al. 2023).

For the implementation of leaf emergence in hydrological models or water balance models, mainly static approaches have been used. These often define annual trends for relevant parameters such as leaf area index (LAI) or vegetated fraction (VCF). These, however, do not allow for variance between years and are therefore not dependent on the actual decisive driving factors. To address this, some common water balance models such as LARSIM (Ludwig & Bremicker 2006) or WaSiM (Schulla 1997) have implemented approaches for the dynamic mapping of phenology, which have hardly been used to date due to inadequate parameterisation approaches. Chen et al. (2023) expanded the SWAT model with a dynamic phenology approach, showing that it improves the model accuracy. However, approaches for the implementation of dynamic phenology in behavioural models are missing but would be necessary to further analyse how hydrological processes and dynamic phenology influence each other.

Although resulting evapotranspiration processes are correct on average when using a static phenology approach, the interannual variability of phenology, i.e. the changing onset of leaf emergence, flowering and senescence, cannot be modelled. A static approach to phenology makes it impossible to correctly simulate the underlying processes, such as evapotranspiration in spring, which depends on leaf emergence (Beaubien & Hamann 2011). Therefore, an implementation of a dynamic phenology approach, i.e. a correct representation of the annual variability of leaf development, would have several advantages for water balance modelling. This could make the mapping of water flux within the system more precise by improving the mapping of transpiration by the plants. This could also improve the correct representation of parameters such as soil moisture and thus have an impact on the modelling of resulting surface runoff. The effects of climate change on leaf emergence and the resulting impact on the water balance in spring could be modelled (Warter et al. 2023). The extension of the vegetation period due to earlier bud burst and later senescence could be depicted. This would allow us to assess their effects on ecological processes. It was also shown that the changing phenology regarding the earlier emergence of leaves greatly influences a catchment's water regime and discharge (Chen et al. 2022, 2023).

Therefore, this study aims to investigate how a dynamic phenology approach can be implemented in a water balance model and whether and how this can improve water balance modelling. For this purpose, we use the dynamic phenology model included in the behavioural water balance model WaSiM. We calibrate it, using a shuffled complex evolution (SCE) algorithm so that the time of leaf emergence in the model comes as close as possible to the time observed in the field. Information on the time of leaf emergence and the evaporation process (MODIS data) is used for this purpose. The optimised model is then validated using Landsat data in relation to the resulting evaporation patterns of the different forest types. For our study, this results in the following objectives and hypotheses: (i) using the data on observed leaf emergence and the evapotranspiration processes according to MODIS, the model-internal dynamic phenology module can be calibrated, whereby (ii) the time of leaf emergence in the model can be dynamically modelled as a function of external parameters and (iii) an improved representation of the evapotranspiration processes for the deciduous and coniferous forests can be achieved.

Stratigraphically, the ridges of the catchment area can be assigned to the Devonian, Lower Devonian, Siegen, Taunus quartzite and Darustwald strata. This means that geologically they consist of quartz sandstone and quartzitic sandstone with intercalations of claystone and siltstone. The valleys are characterised by Quaternary or Pleistocene formations and mainly comprise flowing soils and other sedimentary formations such as slope debris, slope loam and block debris. Cambisols are mainly found on the elevations, and podsolised cambisols and loose brown soils are also found on the two central ridges (Figure 1(b)). The flatter areas between the elevations are flat stagnosols, i.e. areas influenced by stagnant water. Colluvisol-gleysols or pure colluvisols are found along the watercourse with its floodplain and the areas of runoff concentration. Most of the area is covered with forest (Figure 1(c)). The deciduous forest makes up a larger proportion than the coniferous forest. Mixed forest areas only make up a small proportion of the stocking. In the southeastern area, however, there are anthropogenic settlement areas together with agricultural land.

Geological information about the catchment was taken from the ‘Geologische’ ‘Geologische Übersichtskarte 1:300.000 (GÜK300)’ (Landesamt für Geologie und Bergbau 2024), while the soil type information was taken from the ‘Bodenflächendaten im Maßstab 1:50.000 (BFD50)’ (Landesamt für Geologie und Bergbau 2021). Information about the forest types was derived from the European Environment Agency (EEA) (2020), while data for all non-forest land use types were taken from Corine Land Cover (2018). MODIS ETa data (MOD16A2) (Brust et al. 2021) were derived from the NASA Earth data search interface (https://search.earthdata.nasa.gov/search). The mean observation dates for leaf unfolding of Fagus sylvatica in Rhineland-Palatinate, Germany, were taken from the Deutscher Wetterdienst (DWD) (2024). Landsat-8 data (Landsat Collection 2 Level-2 science products) were acquired via the USGS Earth Explorer application (https://earthexplorer.usgs.gov) and included atmospherically corrected surface reflectance (Vermote et al. 2016) and land surface temperature (LST) data (Cook et al. 2014). Gridded station data (Gerlach 2006) for Rhineland-Palatinate were used as meteorological inputs for the WaSiM model and Landsat-8 ETa modelling.

In the model, a rectangular grid of one-dimensional (1D) columns represents the soil. Each column is divided into soil horizons of varying thicknesses, which are also subdivided into several layers each, together with an additional aquifer layer. Van Genuchten parameters and the saturated hydraulic conductivity are used to describe the soil's hydraulic properties. As a result of these varying properties, surface runoff, interflow and groundwater-contributing deep percolation can be generated.

For the model parameterisation, a spatial resolution of 40 m and a temporal resolution of 1 h were chosen as a compromise between reliable model accuracy and computation time. Gridded INTERMET station data (Gerlach 2006) from 01 January 2010 to 31 December 2020 were used as the input time series for meteorological data (temperature, precipitation, radiation, humidity and wind). Following, the preprocessing tool of WaSiM, TANALYS (Schulla 2024b), was used to calculate the required spatial information grids based on the digital elevation model (DEM). These include grids for the slope, exposition, subcatchments, river network, river width and depth. Grids for the aquifer geometries were created, where homogeneous lateral conductivities (⁠ and ⁠) of 5×10⁻⁶ were used. The same value was used for a homogeneous colmation grid. A value of 50 was selected as the threshold for the river network. This threshold value describes how many cells runoff has to be combined in order to form a water body cell in the model. Here, it means that a water body cell is formed after the runoff contribution and concentration from 50 cells. The resulting network showed the best fit with the water body of the catchment. Based on the soil types and land use information, profiles of the individual soils were created. Each profile contains data on thickness, soil type, depth, bulk density, carbonate content, humus content and dry bulk density of the soil's individual horizons. Based on these soil properties, van Genuchten parameters for the model were derived with the pedotransfer function after Wösten et al. (1999). Saturated hydraulic conductivities were taken from the Ad-Hoc AG Boden dataset KA5 (Ad-hoc-AG Boden 2005).

For the phenological simulations, an additional 1D model structure was set up. Instead of simulating the entire catchment, this only simulates a single cell of the catchment grid. This approach was possible because evapotranspiration and phenology essentially only depend on the land use type. The 1D model was then initialised with a podzolic cambisol and an average height of 480 m. One 1D model was set up for the deciduous and coniferous forests.

Satellite-based evapotranspiration data were used as reference data for the model's simulated evapotranspiration. For this, the MODIS ETa data (MOD16A2) (Brust et al. 2021) with a resolution of 500 m for the period from 1 January 2010 to 31 December 2020 was taken. The MODIS data is structured in such a way that each value in the time series represents the sum of the following 8 days. Based on the land use map created for the catchment, five representative MODIS grid cells were selected for the deciduous forest and five representative MODIS grid cells for the coniferous forest, which contain an area of the respective land use type that is as homogeneous as possible. From these five separate time series, a new time series was then created for both land use types, with the mean values of the time series over the respective five grid cells. This meant that a single mean ETa time series was available for each of the two land use types.

The MODIS ETa data were then scaled down by the calculated factor for further simulations.

The mean observation dates for leaf emergence of Fagus sylvatica in Rhineland-Palatinate, Germany (Deutscher Wetterdienst (DWD) 2024), were used as reference data for the exact day of the start of leaf unfolding and are shown in Table 1.

The SCE optimisation algorithm from the SoilHyP package (Andrews et al. 2021) was then used to calibrate the model to the measured emergence dates. Calibration of leaf emergence dates was done for the deciduous forest. The optimisation algorithm was applied via an R script (R Core Team 2023) within RStudio (RStudio Team 2020). For the calibration, the 1D model was initialised with heuristic starting parameter values for ⁠, and ⁠. An initial simulation of the 11-year period from 2010 to 2020 was done. During this, the date of the simulated emergence was written to a separate file for each year. The difference in days between the simulated and the observed emergence day was then calculated for each year. Based on this, the difference sum over the entire 11-year period was determined. This resulted in the cumulative total deviation of the model from the observed values. The value for the total deviation in days was then passed back to the algorithm. The algorithm then generated new parameter values and initialised a new model run with adapted model parameter values. Once a certain number of runs, and therefore sets of parameter values and the associated deviation totals, existed, the algorithm began the optimisation and adjusted the parameter values to minimise the total deviation in days. The result was a set of parameter values that produced the lowest deviation sum for the model and therefore the best fit for leaf emergence estimation.

To check the robustness of the adapted model, a cross-validation of the adjusted model was carried out using a leave-one-out-cross-validation (LOO-CV) approach. From the 11 dates of the calibration period, a new training data set was created that only included 10 of those 11 entries. Based on this data set, the model was trained again using the SCE algorithm approach. With the resulting parameter set, a normal model run was conducted and the deviation in days for the date omitted from the training data set was determined. This was done 11 times so that each date was used once as a validation date. The sum of the deviations of the individual validation dates and the average of these were then calculated. An overview of the individual training sets used, together with their respective calibrated parameters and the day difference values for the validation dates, can be found in Supplementary Table S3.

Once the model had been calibrated to determine the timing of leaf emergence, the ETa curves of the deciduous forest were optimised, using a similar approach. As already shown in Supplementary Table S1, the parameters on which basis the evapotranspiration curve is derived are defined for certain support points (Julian days). For the static parameterisation, parameters are fixed for certain support points and Julian days. However, in the dynamic modelling approach, the timing for the third support point is a fixed Julian day. The timing for the fourth support point is calculated by the dynamic model on the basis of the thermal time. The time for the fifth support point is specified as the difference in days to the start of emergence (in this case the fourth support point). Five parameters were selected for adjustment to optimise the ETa curve for the period around leaf emergence. These are the Julian day for the third support point and the difference between the fourth and fifth points, as well as the LAI values for the third, fourth and fifth support points.

These five parameters were then optimised in a similar way as the parameters for the leaf emergence. The SCE algorithm was used for this approach as well. An initial model run was done with the static model values as starting values for all five parameters, whereby the daily ETa values were written out. These were then converted to 8-day sums in MODIS format. The root mean square error (RMSE) was then determined for the simulated ETa values and the MODIS ETa values over all years. However, only those values were used that were in the period between the 70th day of the year and the 150th day of the year. This was done because the adjustment of the dynamic model has only minimal influence on values outside this interval and the inclusion of identical values from winter periods for the dynamic and static parameterisation outside this range would have influenced the evaluation of the efficiency measure. The calculated RMSE was then passed to the algorithm, which generated new values for the parameters and restarted the model run with these new values. This was repeated until the algorithm converged on the minimised RMSE value. The corresponding parameter set was noted accordingly as the best parameterisation fit.

We decided to not apply this to the coniferous forest as well, as coniferous trees usually lead to significantly less dynamic ETa curves compared with deciduous trees. This is due to the fact that deciduous trees shed all their leaves in fall and form new ones in spring, resulting in a strong contrast in transpiration strength over the year. This is not the case with the conifer species found in the study area, such as Picea abies(L.) H. Karst. It is needled all year round, showing continuous needle senescence of only 1/6 of its needle amount per year (Pretzsch et al. 2012). This is why the ETa curves for the coniferous forest are less dynamic and have a lower amplitude between summer and winter. In addition, due to the smaller difference in foliage cover between summer and winter in coniferous forests, the dynamics in the ETa curves are mainly determined by climatic factors such as temperature and radiation. In deciduous forests, however, a considerable part of the ETa dynamics results from the difference in foliage conditions between summer and winter. As a result, the ETa curves of coniferous forests can only be influenced to a small extent by changing the coniferous leaf growth parameters, whereas in deciduous forests it is possible to exert a much stronger influence via the foliage parameters. We, therefore, used the static parameterisation for the coniferous forest for further analyses.

Landsat-8 evapotranspiration products (50 m) were calculated with the simplified surface energy balance (SSEB) model (Senay et al. 2007, 2013) as described in Casper et al. (2023): In brief, the SSEB model leverages the approximately linear relationship between sensible heat flux and LST within a remote sensing scene to estimate the evaporative fraction (EF) of each image pixel based on the scene's LST distribution. For complex landscapes, the model additionally integrates DEM and vegetation index data (NDVI) into the EF computation to take into account elevation-based LST variation and vegetation cover (Senay et al. 2011). The EF is then pixel-wise multiplied with reference evapotranspiration (ET₀) raster, calculated with standard meteorological station data using the Penman–Monteith equation for a shortgrass crop (FAO-56), to derive pixel-wise ETa estimates.

The calibrated dynamic model was then evaluated to test whether the depiction of the ETa processes could be improved by the adjustment. We decided to use the ratios of the ETa quantities of both land use types, deciduous and coniferous forests, as a test variable and to compare the ratios of the static and dynamic model to observed ratio values. We believe this to be a suitable method for the model evaluation, as the different emergence times of the deciduous forest compared with the static coniferous forest results in a specific ratio between the ETa for the deciduous forest and the coniferous forest for each day. This can therefore be used as an indicator for a correct progression of leaf emergence in the deciduous forest and thus the correctness of the dynamic model.

The ratio of the daily ETa sums of the deciduous forest to that of the coniferous forest was then calculated for each date. This provided one ratio value for the Landsat scenes, the static and the dynamic model for each date. On this basis, the differences between the ratio values for the static model and the Landsat values and between the ratio values for the dynamic model and the Landsat values were determined. This resulted in four deviation values each for the static and the dynamic model. An average value was then formed from these values, which was used to assess the capabilities of the dynamic and static model to depict the correct ETa regimes. In this context, lower deviations indicate better results.

In the first part of the calibration and evaluation of the dynamic model, namely the simulation of the start of leaf emergence, a significant reduction in the deviation from measured dates was achieved (Table 2). The static model achieved an average deviation of 6.09 days, while the dynamic model had an average deviation of 2.36 days. The cumulative deviation over the 11-year period was 67 days for the static model, compared with just 26 days for the dynamic model. There was a pronounced reduction in the deviation for almost all years. Only the years 2011, 2017 and 2020 showed an increase in the deviation in the dynamic model, with the absolute deviation for 2017 being exceptionally large for both models at 11 days for the dynamic and 7 days for the static model.

The cross-validation using the leave-one-out method resulted in a cumulative deviation of 37 days over all 11 years, which corresponds to an average deviation of 3.36 days per year. The day difference values for the cross-validation are similar in range to the calibrated dynamic model (Wilcoxon rank sum test; W = 78, p = 0.253). Despite the cross-validation results in a slight increase in the deviation in days compared with the model calibrated with the full data set, it is still significantly lower compared with the deviation based on the static approach. However, the clear increase in deviation during cross-validation shows that a smaller sample size than the used one would not be sufficient to calibrate the dynamic model correctly.

As part of the adjustment of the support points and the LAI values for the dynamic model, an improvement in both the RMSE and the correlation coefficient, when compared to the MODIS ETa estimates, was achieved for the deciduous forest. The RMSE was decreased from 3.55 for the static model to 3.46 for the dynamic model. The Pearson correlation coefficient was increased from 0.86 to 0.89 as part of the optimisation. The adjusted table with the ETa relevant parameter values of the optimised dynamic model can be found in Supplementary Table S4 for the deciduous forest.

The evaluation of the Landsat scenes and the comparison with the model data for the individual dates showed a significant improvement in the ratios and patterns of the ETa values for three of the four dates in the dynamic modelling approach (Table 3).

The static model shows a 12.1% deviation from Landsat data for the ETa ratio of deciduous to coniferous forests on 27 March 2012, while the dynamic model reduces this deviation to 4.8%, achieving an absolute improvement of 7.33% and a relative improvement of approximately 60%. On 7 April 2018, the static model's deviation is 13.8%, improved by the dynamic model to 4.7%, resulting in a 9.1% absolute reduction and about 66% relative improvement. For 15 April 2015, the static model shows a 20.7% deviation, which the dynamic model reduces to 6.3%, reflecting a 14.4% absolute reduction and a 69.6% relative improvement. However, for 20 April 2017, both models display small deviations of 10.3 and 11.5%, respectively, with no significant change noted.

The determination of the start time of leaf emergence using the dynamic modelling approach in WaSiM shows a significant improvement in the fitting in relation to the observed leaf emergence. Although the approach chosen to determine the emergence time is a relatively simple model based solely on a temperature threshold to be reached and a conceptually implemented dormancy period, the predicted dates match the ones observed in the study area well. This means that this part of the model is generally capable of predicting the time of the start of leaf emergence on the basis of temperature-forcing data.

The analysis revealed that the leaf emergence forecast also deteriorated for a few years, namely for 2011, 2017 and 2020. However, the deviations in 2011 and 2020 are relatively small at 4 and 2 days, respectively, in the dynamic model. Only the year 2017 stands out with a deviation of 11 days. While emergence was not observed in the field until 21 April, the model shows emergence as early as 10 April. One possible explanation for this difference could lie in the way the dynamic approach works within WaSiM. Leaf emergence starts irreversibly as soon as the temperature threshold is reached. It is, therefore, possible that the necessary temperature threshold was reached in the field at this early date, i.e. the model starts leaf emergence at this time, but is then interrupted or delayed by a phase of cold temperatures. For Fagus, for example, it has been shown that cold temperatures and frost can delay the time of leaf emergence (Muffler et al. 2024) and leaf growth after emergence (Rubio-Cuadrado et al. 2021). Leaf loss due to a late frost event after bud break also leads to reduced growth and thus to delayed increase of the evapotranspiration that depends on it (Hufkens et al. 2012; Príncipe et al. 2017). However, the model currently does not reproduce this; once the threshold has been reached, the LAI starts to rise unstoppably according to the parameter curves in the table. An extension of the model would be useful here, which slows down or stops the growth under low temperatures or frost conditions. This would take such temperature drops into account in the simulation of the phenology, as well as the resulting changes in the leaf emergence course.

The improvement in the efficiency measures, as well as in the comparison of the ETa ratios and spatial patterns, shows that the chosen approach is generally suitable for calibrating the dynamic phenology model based on MODIS data. The applied optimisation of the LAI curves and support points resulted in one single specific curve per land use type with the corresponding support points, which was used for the simulation of all years. This is probably the reason for the poor representation of a few years such as 2011 or 2013. If the actual LAI curves in the field are significantly influenced by unusual weather patterns, such as heavy frost events after early budding, and deviate from their typical, average course, then the resulting ETa curves observed in the field will also deviate from the average course. However, since the model uses the same LAI curve every year and does not provide for any modulation of the LAI curves depending on specific environmental parameters, the model cannot depict such deviations from the mean ETa development. Strong interannual variance in the emergence processes of the trees, therefore, leads to reduced applicability of the model. In that case, the resulting simulated ETa curves deviate more strongly from those observed in the field (Wilson & Baldocchi 2000).

An extension of the dynamic approach would make sense here in that no static parameter values are used for the simulation of the ETa, but rather a dynamic determination of these is carried out. Here, a wide variety of approaches could be implemented. Li et al. (2010), for example, have described a dynamic estimation of LAI values based on MODIS data. For an application of that approach, LAI development curves would have to be fitted for each year individually, leading to specific ETa-influencing plant parameter sets for each year, rather than one single parameterisation used for all years. Another approach might also be useful, in which the LAI values are predicted using temperature-dependent regressions for the NDVI for example (Ji & Peters 2004), similar to the temperature sum model for dynamically simulated leaf emergence.

Coupling the WaSiM model with a plant growth model might be a good way as well to achieve dynamic parameterisation of the LAI processes. For example, de Noblet-Ducoudré et al. (2004) used such an approach by coupling the carbon-water-energy model ORCHIDEE (Krinner et al. 2005) with the crop model STICS (Brisson et al. 1998, 2002). As a result, the crop model provided the phenological plant parameters such as the LAI values, which could be used by ORCHIDEE to simulate evapotranspiration, leading to an improved estimation of the evapotranspiration processes. A similar approach was taken by Tsarouchi et al. (2014), who coupled the land surface model JULES (Best et al. 2011; Clark et al. 2011) with the crop growth model InfoCrop to investigate the inter-seasonal land cover changes in evapotranspiration fluxes. Again, evapotranspiration is simulated by the JULES model based on the LAI values contributed by the crop growth model. This shows that such an approach of coupling models that simulate atmospheric and hydrological processes with those that model crop growth is a suitable way to improve the modelling of evapotranspiration.

In the vast majority of existing coupled modelling approaches, however, the plant growth model used is a crop model. This usually is developed for field plants on agricultural land, whereas in our application case, the majority of the areas are not stocked with field plants but with forest stands. Therefore, most of the modelling approaches developed so far would not be applicable to our case. Coupling a water balance model with a forest growth model, therefore, seems to be a sensible approach to dynamically simulate the ETa-relevant plant parameters. Such a forest growth model could be, for example, 4C (‘FORESEE’-FORESt Ecosystems in a Changing Environment), which has been developed for the investigation of long-term changes in the forest structure under changing environmental conditions (Lasch et al. 2005). This would open up the possibility of investigating the effects of a changing precipitation regime on both the water balance and subsequently on the forest itself. This is particularly important because the water balance and the forest influence each other. This means that the simulation of the forest would also benefit from a more dynamic representation of the water balance in terms of its accuracy.

However, the currently implemented features alone open up a wide range of possible applications for the model. For example, it can be used to simulate climate change scenarios by replacing the meteorological inputs with global climate model predictions for different emission scenarios. In this context, the effects of potential changes in phenology behaviour could be analysed by comparing temporal and spatial ETa patterns as well as catchment discharge changes. This is particularly interesting because the module for dynamic phenology is implemented in WaSiM, which facilitates the examination of the effects of changes in the dynamic phenology on the water balance for a given catchment.

With our study, we demonstrate a functional approach to parameterise, calibrate and evaluate the dynamic phenology within a water balance model. We were able to determine the time of the beginning of leaf emergence relatively precisely with our model, reducing the mean deviation from 6.09 to 2.36 days. The calibration of the ETa-determining plant parameters with MODIS ETa time series observations led to an improvement in the fitting of the ETa curves in the deciduous forest. Here, the RMSE was reduced for the deciduous forest from 3.55 to 3.46. The deviation in ETa ratios was improved by 8% on average. This was also confirmed by significantly better fitting ratios of the ETa values of the deciduous forest to those of the coniferous forest. Importantly, the dynamic model showed substantially improved agreement regarding the spatial patterns of ETa when compared against independently estimated Landsat ETa rates. For future work on implementing a dynamic phenology approach, it would be useful not only to dynamically simulate the start of the phenology, but also to dynamically parameterise the course of leaf emergence between the individual years, for example by using time series of remotely sensed vegetation parameters. This could lead to a significant improvement in ETa curves, especially in years that differ more in terms of leaf development behaviour.

We thank Jörg Schulla for his constant support on the handling of the WaSiM model.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project number 426111700 and Forstliche Forschungsförderung Nr. 5.2-04-2023 – Project ‘Klimawald2100 Modul Wald und Wasser’. The publication was also funded/supported by the Open Access Fund of Universität Trier.

M.H. conceptualised the study, did the analysis, wrote and edited the article. C.H. did the remote sensing data analysis, reviewed and edited the article. M.C. conceptualised the study and reviewed the article.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.