An extension of the BROOK90 hydrological model for estimation of subdaily water and energy fluxes

We present an updated version of the BROOK90 hydrological model (B90), which integrates a closed energy and water balance on subdaily time scales. This updated version refines the time-discrete, physically based representations of soil evaporation, interception and transpiration, while also improving the simulation of snowmelt processes on vegetated surfaces. Additionally, the model includes the capability to simulate the sensible heat flux, thereby providing a comprehensive description of the energy balance in conjunction with the water balance for vegetated surfaces.

The model was validated using Integrated Carbon Observation System (ICOS) eddy-covariance measurements from a mature spruce forest at the Anchor Station Tharandt (DE-Tha) in Germany. The subdaily B90 model demonstrates a good agreement with observed 30 min latent and sensible heat fluxes for dry surfaces, while its performance is less accurate for wet surfaces.

Notably, this new version of B90 does not require recalibration; parameter sets from earlier versions remain applicable. The model is well-suited for sampling intervals from 8 h to 1 min, depending on the availability and resolution of the input forcing data. It can be effectively used for various purposes, such as validating flux measurements, gap-filling latent and sensible heat flux data, and performing plausibility checks. However, its primary application is in the study of subdaily water balance dynamics, including processes like dew formation, interception, and fog deposition, which are typically not captured in many other daily-scaled water balance models.

1 Introduction and Motivation

Hydrological modeling plays a crucial role in understanding and managing water resources, predicting hydrological responses to environmental changes, and addressing issues such as flood forecasting and water supply management. The BROOK90 (B90) lumped hydrological model, initially developed in 1978 and subsequently refined by Federer (Federer and Lash, 1978a, b; Federer et al., 2002, 2003), has gained widespread adoption within hydrological research communities due to its robust performance in simulating daily water balance processes. In addition to the original Fortran and Visual Basic implementations (Federer, 2002), the model has seen several adaptations, including forks such as LWF-BROOK90 (Hammel and Kennel, 2001; Schmidt-Walter et al., 2020) and an R-based version (Kronenberg and Oehlschlägel, 2019). A modified Shuttleworth-Wallace evapotranspiration model (Shuttleworth and Wallace, 1985) is incorporated into B90, making it a reliable tool for estimating the water balance of vegetated surfaces. Moreover, this model allows for the separate analysis of transpiration, interception, and evaporation processes at a daily time resolution. Applications of B90 include forest hydrology (Schwärzel et al., 2009), soil monitoring systems (Hohenbrink et al., 2024; Luong et al., 2023; Vorobevskii et al., 2024), and even flood forecasting (Luong et al., 2021). Although B90 was originally designed for small plot-scale simulations, it has been successfully applied to catchments exceeding 100 km² (Ulker and Buyukyildiz, 2023; Vorobevskii et al., 2020).

Despite its success in accurately simulating atmosphere-plant-soil water fluxes at the point scale, the original B90 model operates at a daily output time step and does not close the energy balance, focusing solely on the water mass balance. This limitation restricts its ability to capture subdaily dynamics of hydrological processes. While Federer (2002) asserted that subdaily output time-steps would not be suitable for B90, we aim to demonstrate that this is not the case. Both mass and energy exchanges between the land surface and the atmosphere are critical processes that must be addressed for accurate water balance estimations in hydrological models operating on subdaily scales.

Recent advancements in meteorological forcing data, such as high-resolution micro-meteorological measurements (Pastorello et al., 2020) and reanalysis datasets like ERA5-Land (Muñoz-Sabater et al., 2021), offer new opportunities to incorporate additional features of hydro-meteorological processes into water balance estimation. For example, rainfall partitioning can be significantly influenced by antecedent soil moisture, which exhibits substantial temporal variation even on a subdaily scale (Taylor et al., 2012). Additionally, fog and dew processes, which can contribute notably to water inputs in specific ecosystems, are often neglected in daily models (Parlange et al., 1995; Wilson et al., 2000). Dew formation is a nocturnal process driven by radiative cooling, and neglecting the associated energy loss during the night can lead to underestimations of dew formation (Monteith, 1957; Jacobs et al., 1999, Körner et al., 2020). Similarly, fog deposition requires precise energy balance considerations to be accurately modeled (Calder, 1996; Gash et al., 1980). Furthermore, event-driven processes such as rain and snow interception, which affect both water and energy fluxes, necessitate a model capable of operating at subdaily time resolutions.

To overcome these limitations, we present an enhanced version of the BROOK90 model, incorporating a new algorithm that ensures closure of both the energy and water balances, enabling simulations and output at subdaily time scales. This development builds upon the established Shuttleworth-Wallace evaporation model and extends its application to subdaily processes. By partitioning evapotranspiration into transpiration, evaporation, and interception processes, the updated B90 model offers a more comprehensive and accurate representation of hydrological dynamics. This advancement is particularly important for applications in ecosystem management, agricultural planning, and water resource management, where understanding the detailed subdaily dynamics of water and energy fluxes is crucial.

2 Location and eddy-covariance measurement data

Data from the Integrated Carbon Observation System (ICOS) Anchor Station Tharandt (DE-Tha) station in Saxony, Germany, collected during the 2023-2024 period at a 30 min temporal resolution (Bernhofer et al., 2024) were selected to demonstrate the performance and evaluate the subdaily B90 model (Fig. 1). This ICOS station's data has been extensively used in numerous studies and is recognized for its high quality (Moderow et al., 2021). The site is dominated by a 140-year-old spruce stand with an average canopy height of 30 m. The tower's footprint encompasses 87 % coniferous forest (72 % Picea abies) and 13 % deciduous forest (10 % Larix decidua). The forest has been under continuous management since 1811, involving planting, thinning, and liming practices. The typical root depth of the spruce trees ranges from 30 to 40 cm.

Figure 1 Location of the ICOS DE-Tha station: (a) overview map, (b) satellite image (Bing Satellite © Microsoft, 2024) of the site with meteorological station and eddy-covariance tower, (c) photo of the eddy-covariance tower above the forest stand (© Chair of Meteorology, TU Dresden).

The site is characterized by silty Podzol soils with a relatively high stone content (10 %–20 %), derived from periglacial deposits of rhyolitic debris and loess, resulting in a highly heterogeneous composition. The eddy-covariance measurement system is installed atop a 42 m tower, equipped with ultrasonic anemometers and closed-enclosure gas analyzers, providing data at a frequency of 25 Hz. In addition, the site is outfitted with instruments to measure radiation components, air temperature and humidity, precipitation, and to conduct soil and biomass observations.

For the B90 model forcing, the following measurements were used as input: minimum and maximum air temperature, shortwave incoming radiation, precipitation, wind speed, relative humidity, and ground heat flux. Model validation was carried out using measured latent heat flux (λE) and sensible heat flux (H), with ICOS quality flag zero. The B90 model parameters for vegetation and soil were derived from Vorobevskii et al. (2022). The year 2023 was used as the spin-up period for the model and was not included in the evaluation phase.

3 Extension of the model: energy balance components and closure

Evapotranspiration in the B90 model is based on an enhanced version of the Shuttleworth-Wallace approach (Shuttleworth and Wallace, 1985). It is a dual-source model (Fig. 2), which explicitly differentiates between wet and dry surface conditions. Derived from the Penman-Monteith model (Allen et al., 1998), the Shuttleworth-Wallace model retains the three key equations of the Penman-Monteith framework – namely, the turbulent fluxes of latent heat (λE) and sensible heat (H), as well as the surface temperature (T_o). In this section, we present the newly implemented equations in subdaily B90 for the estimation of sensible heat flux (H).

Figure 2 Scheme of the minimal sensible heat flux adapted after Shuttleworth and Gurney (1990) in subdaily BROOK90.

The minimal sensible heat flux H_min in the subdaily B90 model, as depicted in Fig. 2, represents the sum of five processes originating from two sources: the canopy and the soil. The canopy can be either dry or wetted by precipitation, but a distinction between fluxes from liquid and solid precipitation is not made for the canopy. Depending on the occurrence of a precipitation event or residual water from a prior event (i.e., interception), the leaf surface may be fully or partially covered by water.

Similarly, the soil can be covered by either rain or snow during a precipitation event. For the soil, the differentiation between rain and snow is critical, as snow forms its own storage when accumulated, whereas rain either infiltrates into the soil or runs off. Each of these five flux components contributes to the overall heat flux at screen height. Consequently, the minimal sensible heat flux H_min is estimated according to Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{H}_{\min }& =\left(\mathrm{1}-w\right)\cdot H_{\mathrm{dry}}^{\mathrm{Canopy}}+w\cdot H_{\mathrm{wet}}^{\mathrm{Canopy}}+H^{\mathrm{Snow}}\\ & +\left(\mathrm{1}-s\right)\cdot \left(\left(\mathrm{1}-w\right)\cdot H_{\mathrm{dry}}^{\mathrm{Soil}}+w\cdot H_{\mathrm{wet}}^{\mathrm{Soil}}\right)\end{array}\end{array}$$

where w is the portion of wetted surface in the canopy and on the ground. It is defined as the quotient between the actual intercepted rain (IRVP) and snow (ISVP) divided by the potential interception rate of the canopy (PINT). In case the quotient gets larger than one, it is reduced like shown in Eq. (2):

$$\begin{array}{}\text{(2)}& w=\min \left({\displaystyle \frac{\mathrm{IRVP}+\mathrm{ISVP}}{\mathrm{PINT}}},\mathrm{1}\right)\end{array}$$

Equation (1) possess an additional trigger for snow s. In case of snow on the soil surface evaporation of this source is estimated and the snow on the ground trigger s is set to one, otherwise it is zero:

$$\begin{array}{}\text{(3)}& s=\left\{\begin{array}{c}\mathrm{0},(\mathit{\lambda }{E}^{\mathrm{Snow}}\le \mathrm{0})\\ \mathrm{1},(\mathit{\lambda }{E}^{\mathrm{Snow}}>\mathrm{0})\end{array}\right.\end{array}$$

The sensible heat flux from snow pack on the ground is estimated according to:

$$\begin{array}{}\text{(4)}& H^{\mathrm{Snow}}=-k\cdot {\displaystyle \frac{\mathit{\rho }}{\mathit{\gamma }}}\cdot {\displaystyle \frac{{\mathrm{VPD}}^{\mathrm{Snow}}}{r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{s}}}}\end{array}$$

where VPD^Snow is the vapour pressure deficit above snow pack on the ground, ρ is the air density, γ is the psychrometer constant, $r_{\mathrm{a}}^{\mathrm{a}}$ is the atmosphere aerodynamic resistance, $r_{\mathrm{a}}^{\mathrm{s}}$ ground aerodynamic resistance, k is a multiplier to fix overestimation of snow evaporation introduced by Federer (2002).

The sensible heat flux from dry canopy is defined as:

$$\begin{array}{}\text{(5)}& H_{\mathrm{dry}}^{\mathrm{Canopy}}={\displaystyle \frac{\left(\mathrm{AE}-{\mathrm{AE}}_{\mathrm{SUBS}}\right)\mathit{\gamma }\left(r_{\mathrm{s}}^{\mathrm{c}}+r_{\mathrm{a}}^{\mathrm{c}}\right)-\mathrm{VPD}\mathit{\rho }}{r_{\mathrm{a}}^{\mathrm{c}}\left(\mathit{\gamma }+\mathrm{\Delta }\right)+\mathit{\gamma }{r}_{\mathrm{s}}^{\mathrm{c}}}}\end{array}$$

where AE is the available energy at screen height, AE_SUBS is the part of energy reaching the ground, $r_{\mathrm{a}}^{\mathrm{c}}$ is the canopy aerodynamic resistance, $r_{\mathrm{s}}^{\mathrm{c}}$ is the canopy surface resistance after Shuttleworth and Gurney (1990) and Stewart (1988), Δ is the rate of change of saturation specific humidity with air temperature.

The sensible heat flux from wet canopy is defined as:

$$\begin{array}{}\text{(6)}& H_{\mathrm{wet}}^{\mathrm{Canopy}}={\displaystyle \frac{\left(\mathrm{AE}-{\mathrm{AE}}_{\mathrm{SUBS}}\right)\mathit{\gamma }\left(r_{\mathrm{a}}^{\mathrm{c}}\right)-\mathrm{VPD}\mathit{\rho }}{r_{\mathrm{a}}^{\mathrm{c}}\left(\mathit{\gamma }+\mathrm{\Delta }\right)}}\end{array}$$

The sensible heat flux from dry soil is defined as:

$$\begin{array}{}\text{(7)}& H_{\mathrm{dry}}^{\mathrm{Soil}}={\displaystyle \frac{\left({\mathrm{AE}}_{\mathrm{SUBS}}\right)\mathit{\gamma }\left(r_{\mathrm{s}}^{\mathrm{s}}+r_{\mathrm{a}}^{\mathrm{s}}\right)-\mathrm{VPD}\mathit{\rho }}{r_{\mathrm{a}}^{\mathrm{s}}\left(\mathit{\gamma }+\mathrm{\Delta }\right)+\mathit{\gamma }{r}_{\mathrm{s}}^{\mathrm{s}}}}\end{array}$$

where $r_{\mathrm{s}}^{\mathrm{s}}$ is the soil surface resistance after Shuttleworth and Wallace (1985).

The sensible heat flux from wet soil is defined as:

$$\begin{array}{}\text{(8)}& H_{\mathrm{wet}}^{\mathrm{Soil}}={\displaystyle \frac{\left({\mathrm{AE}}_{\mathrm{SUBS}}\right)\mathit{\gamma }\left(r_{\mathrm{a}}^{\mathrm{s}}\right)-\mathrm{VPD}\mathit{\rho }}{r_{\mathrm{a}}^{\mathrm{s}}\left(\mathit{\gamma }+\mathrm{\Delta }\right)}}\end{array}$$

It is important to note that the calculation of the canopy surface resistance $r_{\mathrm{s}}^{\mathrm{c}}$ in B90 differs from the Jarvis (1976) model in that Federer (2002) defined the resulting transpiration term in B90 as “potential” transpiration, neglecting the dependency on soil water availability and CO₂ concentration. As a result, B90 simulates the actual water supply rate, which is controlled by the water potential gradient and plant resistance.

It could be noticed that we have so far formulated only a minimal sensible heat flux H_min. This requires some explanation. B90 was originally developed to represent the water balance, focusing primarily on the mass balance. Consequently, certain components of input energy were not explicitly considered, as energy components were not intended to be modeled outputs. The energy balance is not explicitly closed in the original B90 model. This is because B90, following the Shuttleworth-Wallace approach, estimates a potential latent heat flux, which is limited by the actual available water. Therefore, any portion of the potential latent heat flux that cannot be utilized disappears. This disappearance needs to be addressed in order to close the energy balance in B90. Since it is assumed that this energy is converted into sensible heat, we refer to this process as the redistribution of energy within the model.

We identified four processes in B90 that require such redistribution to conserve the energy balance. At each simulation time step, this redistribution portion ΔH begins with zero value. This occurs for two reasons: first, energy storage is not considered in the model, and second, not all of the identified processes contribute at every time step. The following four steps summarize the amount of redistributed energy ΔH for a single modeling time step:

• Step 1. Redistributed energy from transpiration.

  The main processes in B90 where energy might “disappear” is transpiration, where a “potential” transpiration rate is calculated but the actual amount is limited by the water availability in the soil. At each time step their difference is added to ΔH:

  $$\begin{array}{}\text{(9)}& \mathrm{\Delta }H=\mathrm{\Delta }H+\mathit{\lambda }{E}_{\mathrm{potential}}^{\mathrm{Canopy}}-\mathit{\lambda }{E}_{\mathrm{actual}}^{\mathrm{Canopy}}\end{array}$$

  with $\mathit{\lambda }{E}_{\mathrm{potential}}^{\mathrm{Canopy}}$ estimated after Shuttleworth-Wallace and $\mathit{\lambda }{E}_{\mathrm{actual}}^{\mathrm{Canopy}}$ as actual water flux limited by the available water.

• Step 2. Redistributed energy from snow pack on the ground.

  The second process is evaporation from snow, where a “potential” evaporation rate is calculated, but the actual water vapour uptake is limited by the remaining water after melted runoff. Therefore, remaining parts of latent energy, which are not assigned to the latent heat flux need to be attributed to the sensible heat flux from the snowpack.

  $$\begin{array}{}\text{(10)}& \mathrm{\Delta }H=\mathrm{\Delta }H+\mathit{\lambda }{E}_{\mathrm{potential}}^{\mathrm{Snow}}-\mathit{\lambda }{E}_{\mathrm{actual}}^{\mathrm{Snow}}\end{array}$$

  with $\mathit{\lambda }{E}_{\mathrm{potential}}^{\mathrm{Snow}}$ as latent heat flux equivalent to Eq. (4) after Federer (2002) and $\mathit{\lambda }{E}_{\mathrm{actual}}^{\mathrm{Snow}}$ as actual water flux limited by the available water.

• Step 3. Redistributed energy from snow covered soil.

  In case of a snow cover on the soil the estimated evaporation rates from dry and wet soil surfaces, as well as the sensible heat fluxes are set to zero, although energy is assigned to these processes at each time step. This energy is distributed into the overall sensible heat flux by:

  $$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{\Delta }H& =\mathrm{\Delta }H+\left(\mathrm{1}-w\right)\cdot H_{\mathrm{dry}}^{\mathrm{Soil}}+w\cdot H_{\mathrm{wet}}^{\mathrm{Soil}}\\ & +\left(\mathrm{1}-w\right)\cdot \mathit{\lambda }{E}_{\mathrm{dry}}^{\mathrm{Soil}}+w\cdot \mathit{\lambda }{E}_{\mathrm{wet}}^{\mathrm{Soil}}\end{array}\end{array}$$

  with dry and wet soil evaporation $\mathit{\lambda }{E}_{\mathrm{dry}}^{\mathrm{Soil}}$, $\mathit{\lambda }{E}_{\mathrm{wet}}^{\mathrm{Soil}}$ estimated after Shuttleworth-Wallace.

• Step 4. Redistributed energy from dry soils on wetted by precipitation event.

  In case ground evaporation is limited by the available soil water content and a precipitation event occurs at a simulation time step. The energy difference between potential soil evaporation and actual soil evaporation is redistributed to the sensible heat flux.

  $$\begin{array}{}\text{(12)}& \mathrm{\Delta }H=\mathrm{\Delta }H+\left(\mathrm{1}-w\right)\cdot \left(\mathit{\lambda }{E}_{\mathrm{dry},\mathrm{potential}}^{\mathrm{Soil}}-\mathit{\lambda }{E}_{\mathrm{dry},\mathrm{actual}}^{\mathrm{Soil}}\right)\end{array}$$

  with dry soil evaporation $\mathit{\lambda }{E}_{\mathrm{dry},\mathrm{potential}}^{\mathrm{Soil}}$ after Shuttleworth-Wallace and $\mathit{\lambda }{E}_{\mathrm{dry},\mathrm{actual}}^{\mathrm{Soil}}$ as actual water flux limited by the available water.

  There seems to be no general rule, which process contributes the largest amount of energy for redistribution. In our case study the transpiration process had the largest amount, but at other less vegetated sites or in colder climates also snow evaporation might possess more importance.

  Based on the aforementioned equations the actual sensible heat flux H at a simulation time step is defined as:

  $$\begin{array}{}\text{(13)}& H=H_{\min }+\mathrm{\Delta }H\end{array}$$

  By this redistribution approach the energy balance is closed in subdaily B90. Fortunately, all presented model extensions do not need any new parameters or their adjustments. Therefore, all parameters of the original B90 can be used and should be set according to natural site conditions.

For the technical implementation of the model upgrades mentioned above we took R-Version of the B90 model and upgraded it to a new “subdaily version” (Kronenberg and Vorobevskii 2025).

4 Results and discussion

4.1 Two examples of subdaily water and energy balance component distributions after rainfall events

The updated version of the B90 model allows for the simulation of all water balance components, particularly evapotranspiration, under varying conditions (wet and dry) at subdaily time steps. All output variables previously available in the original B90 model are now also simulated at subdaily intervals. Additionally, new variables have been introduced to represent energy and turbulent fluxes. Here and throughout the manuscript, subdaily time steps refer to 30 min intervals, corresponding to the temporal resolution of the available observational data. As demonstrated in Fig. 3, the model estimates transpiration, evaporation from the ground, and evaporation from the interception storage. Moreover, the model is capable of simulating condensation on vegetative surfaces, along with the subsequent evaporation of intercepted condensed water.

Figure 3 Sample time-series of evaporation components estimated by subdaily B90 and soil water content (simulated and measured): (a) and (c) 29 June–2 July 2024 (9.6 mm on total), (b) and (d) 4–7 December 2024 (22.5 mm in total). Both examples are no snow events.

Figure 3a illustrates this process on the morning of 1 July, where the evaporation rate from intercepted rain becomes slightly negative. A more pronounced example of condensation is observed in Fig. 3b, around midnight on the 5 and 7 December. However, these examples raise additional points that warrant further discussion.

Figure 3c and d illustrate the theta anomaly, defined as the difference between the average θ values of the two or four preceding daily periods and the actual θ values observed in situ and estimated by the B90 model at different soil depths. The theta anomaly is used as a normalized measure to highlight soil water response to precipitation events, as absolute values may vary significantly across different conditions.

In Fig. 3c, the precipitation events of 9.6 mm in total on 30 June–2 July 2024 are shown, with the model capturing a clear response in the soil matrix. In contrast, the in situ measurements at various depths show almost no discernible reaction to the event. The model simulates a drying process for the pre-event soil moisture, yet this is not reflected in the observed data. Several factors may explain this discrepancy: first, an actual interception storage capacity could be larger than parameterized in the model, although this is unlikely due to relatively high sum of precipitation event in a given time interval; second, the presence of an unaccounted-for humus layer above the soil, which can possess a significant storage capacity (Floriancic et al., 2023); or third, the possibility that soil moisture sensors, particularly under dry conditions, exhibit delayed responses or fail to detect moisture changes due to preferential flow bypassing the sensors. These sensor-related issues are further highlighted in Fig. 3d for December with 22.5 mm rain event, where the soil's response to precipitation is shown to produce similar peaks at various depths, consistent with the model's simulation, indicating that the observed anomalies may be attributed to sensor limitations under specific conditions.

Additionally, the latent and sensible heat fluxes are now explicitly simulated in the subdaily B90 model. Exemplary results for the previously discussed events are shown in Fig. 4. A notable observation is that the modelled fluxes exhibit less temporal variability compared to the flux measurements. This is expected, as the observed fluxes are inherently more variable due to the high-frequency nature of the measurements, while the model outputs are influenced by less frequently measured meteorological variables. Furthermore, the model tends to underestimate low or negative fluxes, which can be identified as a general trend. This underestimation is primarily attributed to the systematic underestimation of negative net radiation (Rn) within the B90 model. An improvement in the estimation of longwave radiation within the model is anticipated to address this issue. On the positive side, the introduction of subdaily processes in the model allows for the simulation of condensation events, as well as the representation of oasis effects – where advective warm and dry air leads to negative sensible heat fluxes under specific atmospheric conditions. These effects are visible in Fig. 4b and d.

Figure 4 Sample time-series of energy balance components measured at the ICOS DE-THA site (a and c) and estimated by subdaily B90 (b and d) for 29 June–2 July and 4–7 December 2024.

Table 1 summarizes the energy flux components for June and December 2024. The differences between observed and modeled turbulent fluxes are relatively small, indicating good agreement between the two datasets. In the B90 model, the change in energy storage is assumed to be zero, as no explicit heat storage is implemented. Consequently, the change in storage (ds) is zero at every simulation time step. In December, an underestimation of negative net radiation (Rn) is observed, with both the observed Rn −G and sensible heat flux (H) being lower than the model estimates. However, there is also an unexplained discrepancy, where the observed Rn −G is consistently lower than the modeled values. This issue will be addressed in further discussion.

Table 1 Mean 30 min energy flux components for June and December 2024.

4.2 Validation of energy fluxes with ICOS measurements for 2024

The simulated fluxes were validated against ICOS observations flagged with the highest data quality (quality flag 0 – measured directly). Figure 5 displays all 30 min fluxes, differentiated by vegetation period (VP – April to October) and winter period (WP – November to March). Among the components, latent heat flux (λE) exhibits the weakest performance. This reflects both the inherent complexity of turbulent fluxes and the associated challenges in their accurate modelling and measurement. The observed λE data also contain notable variability, which may be partially attributable to measurement noise. Linear regression analysis for the vegetation period reveals that λE is underestimated by approximately 36 % in the subdaily B90 model, with an acceptable coefficient of determination (R² = 0.59). Model performance further deteriorates during the winter period, characterized by weaker fluxes and similarly pronounced underestimation. Several factors likely contribute to the observed spread in λE:

• the complexity of the underlying evapotranspiration processes,

• the simplified nature of the B90 model, which relies on parameterized rather than turbulence-resolving approaches,

• limitations of eddy covariance measurements under stable atmospheric conditions, and

• uncertainties associated with the site-specific parameterization.

Model performance for λE could potentially be improved through a subdaily calibration of vegetation-related parameters tailored to the study site. In contrast, the modelled sensible heat flux (H) performs much better, showing a 32 % underestimation but a higher correlation with observations (R² = 0.76) over 30 min intervals. This improved agreement likely reflects the relatively simpler dynamics governing H, which depend primarily on air temperature gradients and wind, in contrast to λE, which is additionally constrained by water availability and governed by energy partitioning processes within the model. When H and λE are plotted against the available measured energy (Rn −G), results are satisfactory, showing minimal energy underestimation and a high degree of explained variance. Moreover, the combined turbulent fluxes (H+ λE) from the B90 model exhibit good agreement with the corresponding ICOS observations, further supporting the model's capability to simulate total energy fluxes with reasonable accuracy.

Figure 5 Evaluation of 30 min-resolution energy fluxes from B90 model for year 2024 with eddy-covariance measurements (with quality flag “0” – observed): (a) latent heat, (b) sensible heat, (c) energy balance closure, (d) total energy flux.

Figure 6 and Table 2 presents the cumulative sensible (H) and latent heat fluxes (λE) simulated by the B90 model in comparison with the interquartile range (25th to 75th percentiles) of the processed ICOS observations. The modelled H aligns closely with the median (50th percentile) of the ICOS data, whereas lies within the interquartile range, specifically between the 25th and 50th percentiles of the observed data, indicating a systematic underestimation potentially attributable to the selected site-specific parameterization. The relatively large interquartile range (IQR) in Fig. 6 is expected and reflects the inherent variability of eddy-covariance measurements. High-frequency observations (25 Hz) show considerable variance due to strong winds at 40 m height and the large, heterogeneous tower footprint (up to ∼ 300 m). When fluxes are aggregated to 30 min intervals and cumulative values are calculated, this variability propagates over time, resulting in the observed wide IQR. During the initial quarter of the vegetation period, the cumulative energy fluxes for both H and λE show similar trends between model and observations. However, from that point onward, the model begins to increasingly underestimate λE, a trend that continues until the end of the year. In contrast, the underestimation of H diminishes with the onset of winter, likely due to the reduced role of vegetation-mediated processes in cooler months. The most notable discrepancy is the persistent underestimation of net radiation (Rn) in the ICOS dataset relative to the cumulative turbulent fluxes (Fig. 6c). We could not identify a clear explanation for this mismatch. Notably, the cumulative sum of available energy (Rn − G) from the ICOS data is approximately 18 % lower than both the observed and modelled cumulative turbulent fluxes. This discrepancy may reflect unaccounted energy storage components (keeping in mind that the time-series are discontinuous – only values with best quality flag are filtered, so we don't know if and when the storage gets emptied) within the canopy or soil, which are not explicitly considered in either the measurement gap-filling or the current model configuration.

Figure 6 Cumulated observed and modelled 30 min energy fluxes for the year 2024 with quality flag “0” (observed): (a) latent heat, (b) sensible heat, (c) total available energy.

Table 2 Subdaily quantile errors between observed and simulated latent and sensible heat fluxes by month (2024) with quality flag “0” (observed).

Figure 7 Comparison of different weather conditions for 30 min-resolution latent heat flux: (a) dry weather (transpiration and soil/snow evaporation), (b) wet and snowy weather (snow interception), (c) wet and rainy weather (rain interception).

Figure 8 Evaluation of daily energy fluxes from B90 model for the for year 2024 with eddy-covariance measurements: (a) latent heat, (b) sensible heat, (c) total energy flux.

In the B90 model framework, dry conditions are defined as time steps without precipitation, during which both the canopy and ground surfaces are free of water. In contrast, wet conditions refer to periods when water or snow is stored on the canopy or ground surfaces, but no precipitation is actively occurring. Our results indicate a clear difference in the model's performance in simulating latent heat flux (λE) under dry versus wet conditions. As illustrated in Fig. 7, model performance is highest under dry conditions. During these periods, λE is primarily driven by transpiration, suggesting that B90 effectively captures this component of the evapotranspiration process. The wider scatter observed in the data may be attributed to several factors, including the omission of CO₂ effects on stomatal regulation, oversimplified or inaccurate temporal parameterization of the leaf area index (LAI), and the lack of representation of water storage in plant stems. In contrast, model performance under wet conditions – such as those associated with rain or snow – shows little to no correlation with observed latent heat fluxes from ICOS eddy covariance measurements. Precipitation type (rain or snow) was distinguished using an air temperature threshold. The poor agreement under these conditions may be due to limitations in the interception module of B90 and/or reduced measurement accuracy of eddy covariance systems during precipitation events. It is well documented that water-covered anemometer sensors and flux towers can introduce biases under such conditions. Further investigation into interception and throughfall processes, supported by targeted field observations, could help to elucidate these discrepancies and improve model representation.

Figure 9 Comparison of evaporation components between original daily and aggregated subdaily B90 models for the year 2024: (a) transpiration, (b) snow and rain interception, (c) evaporation from soil and snow.

The aggregation of 30 min simulation outputs to daily values yields highly satisfactory agreement with observations, as illustrated in Fig. 8. Among the energy flux components, latent heat flux (λE) exhibits the lowest performance, with a coefficient of determination (R²) of 0.79. In contrast, sensible heat flux (H) and the combined flux (λE + H) show stronger correlations with measured data, indicating a robust representation of the overall energy exchange. The systematic underestimation of simulated energy fluxes – particularly λE – may be attributable to the omission of energy storage terms within the canopy. Incorporating such storage processes in future model developments could enhance the accuracy of energy balance closure.

Additionally, we compared the 24 h aggregated latent heat flux (λE) values from the subdaily B90 simulations with the corresponding daily λE outputs from the original B90 model (Fig. 9). As shown in Fig. 9, the same parameter set was applied to both the original and subdaily B90 model configurations. This consistency indicates parameter stability across varying simulation time steps. Additional tests were conducted using a range of temporal resolutions (0.5, 1, 2, 3, 4, 6, 8, and 12 h). These tests revealed no significant improvement in model performance associated with a specific aggregation interval. However, model accuracy begins to degrade when subdaily B90 is applied with time steps exceeding 8 h, particularly when compared to the original daily B90 results. This decline is attributed to the inherent design of the original B90, which incorporates a fixed diurnal (day–night) partitioning of mass fluxes. Consequently, subdaily B90 simulations using a 12 h interval (e.g., split at midnight and midday) cannot adequately resolve diurnal dynamics and fail to capture the distinct day and night flux components.

Therefore, we recommend the application of the subdaily B90 model for simulation intervals of 8 h or less. This recommendation is based on the structure of the original model, which incorporates a day–night separation in the evapotranspiration process, resulting in improved performance at coarser temporal resolutions – particularly at daily (24 h) time steps. However, this does not imply that daily aggregates derived from subdaily B90 simulations are inferior. On the contrary, the model preserves the integrity of process-based simulations across temporal scales, ensuring that high-resolution outputs can be reliably aggregated to longer time steps without loss of data quality.

5 Conclusions

We presented an enhanced version of the BROOK90 hydrological model capable of simulating water and energy fluxes at subdaily temporal resolution for vegetated land surfaces. In addition to implementing subdaily time steps, we explicitly close the energy balance within the new B90 framework. Model validation for the selected study site demonstrates good agreement between eddy covariance measurements and simulated fluxes at 30 min intervals.

While the updated model shows promising results, several simplifications were necessary to avoid introducing additional parameters. Notably, processes such as energy and water storage within the canopy and soil, which are relevant at subhourly scales (i.e., Federer, 2002), are currently omitted. Their inclusion would require new parameters that have not yet been incorporated in this version.

Moreover, as turbulent exchange becomes increasingly important with finer temporal resolutions, it is essential to account for atmospheric stability in estimating aerodynamic resistances. This aspect is especially relevant for correctly simulating turbulent fluxes (Banerjee et al., 2017).

Despite these limitations, the subdaily implementation of B90 provides valuable opportunities to investigate long-term and climatological variations in water and energy balances at finer temporal scales. Unlike most hydrological models that operate on a daily basis, subdaily B90 enables the detailed study of high-frequency processes such as dew formation, interception, and fog deposition.

Appendix A: Description of the instruments and measured variables at the ICOS DE-Tha station used in this study

• Eddy-covariance system – GA_CP-LI-COR LI-7200, SA-Gill,

• HS-50, Grill R3-50, 42 m height,

• Air temperature and relative humidity – RHTEMP-Vaisala HMP45, 40 m height,

• Precipitation – Pluvio2 (OTT Hydromet), 1 m height,

• Wind speed - SA-Gill HS-50, Grill R3-50, 42 m height,

• Shortwave radiation – CNF4 (Kipp&Zonen), 37 m height,

• Soil heat flux – SOIL_H-Hukseflux HFP01SC, 5 cm depth,

• Soil moisture – SWCTEMP-Campbell CS65X, 2–50 cm depth.

Appendix B: Additional validation of the subdaily BROOK90 model using eddy covariance data from other stations

Here we present validation results of the subdaily BROOK90 model using additional ICOS-class eddy-covariance stations located in the vicinity of the DE-Tha station and representing different land-cover types. Table B1 summarizes the location, land cover, and available time period of each station used for the additional validation. The model parameterisation for these stations was adopted from Vorobevskii et al. (2022). The evaluation was based on the comparison of simulated and observed latent heat flux (λE) and sensible heat flux (H) for 30 min time-scale using the Kling–Gupta efficiency (KGE) metric and its components. A half-year model warm-up period was applied.

Table B2 presents the validation results using all available eddy-covariance observations. This comparison provides an overall assessment of model performance across the full available time series, including measured and gap-filled or lower-quality observations.

Table B3 presents the same validation metrics, but only for observations with quality flag 0, representing directly measured high-quality flux data. This stricter filtering allows the model performance to be assessed against the most reliable observational subset.

Table B1 Summary of the ICOS eddy-covariance stations used for additional validation.

Table B2 Validation results – all available data.

Table B3 Validation results – only data with quality flag “0”.