From weather data to river runoff: using spatiotemporal convolutional networks for discharge forecasting

The quality of river runoff determines the quality of regional climate projections for coastal oceans or other estuaries. This study presents a novel approach to river runoff forecasting using convolutional long short-term memory (ConvLSTM) networks. Our method accurately predicts daily runoff for 97 rivers within the Baltic Sea catchment by modeling runoff as a spatiotemporal sequence defined by atmospheric forcing. The ConvLSTM model predicts river runoff with an accuracy of ±5 % when compared to the hydrological model. Compared to more complex process-based hydrological models, ConvLSTM networks offer fast processing times and easy integration into climate models, demonstrating their potential as a powerful tool for climate simulation and water resource management.

1 Introduction

River runoff is a key component of the global water cycle as it comprises about one-third of the precipitation over land areas (Hagemann et al., 2020), making accurate runoff forecasting essential for effective water resource management, particularly over extended periods (Fang et al., 2019; Tan et al., 2018). In addition to short-term forecasting, long-term projections of river runoff are vital for climate change studies, projecting flooding and droughts over global and river basins (Cook et al., 2020). These studies calculate river runoff using a land model incorporating a hydrological model within a coupled Earth system model (ESM) (Wang et al., 2022). In the absence of a fully coupled ESM, river runoff as input for ocean models can be created using hydrological models such as Hydrological Discharge (HD) (Hagemann et al., 2020) or HYdrological Predictions for the Environment (HYPE) (Lindström et al., 2010). Hydrological models represent a process-based approach where the water balance is calculated using hydrological processes (e.g., snow, glaciers, soil moisture, or groundwater contribution). These models are complex forecasting tools that are widely utilized, such as high-resolution multi-basin models applied across Europe (Hundecha et al., 2016).

The second approach to projecting river runoff employs data-driven models, such as when calculating river runoff as the difference between precipitation and evaporation over a catchment area, with an integrated statistical correction (Meier et al., 2012). With the recent rise in machine learning in climate research, various data-driven model architectures have been explored for river runoff forecasting. Common approaches include feed-forward artificial neural networks; support vector machines; adaptive neuro-fuzzy inference systems; and, notably, long short-term memory (LSTM) neural networks. LSTM networks have gained traction for long-term hydrological forecasting due to their excellent performance (Humphrey et al., 2016; Huang et al., 2014; Ashrafi et al., 2017; Liu et al., 2020; Fang and Shao, 2022; Kratzert et al., 2018). LSTM networks, first introduced by Hochreiter and Schmidhuber (1997), are an evolution of classical recurrent neural networks (Sherstinsky, 2020). A significant advantage of an LSTM network's architecture is the memory cell's ability to retain gradients. This mechanism addresses the vanishing gradient problem, where, as input sequences elongate, the influence of initial stages becomes harder to capture, causing gradients of early input points to approach zero. LSTM networks have shown stability and efficacy in sequence-to-sequence predictions. However, a limitation of LSTM networks is their inability to effectively capture two-dimensional structures, which is an area where convolutional neural networks (CNNs) excel (Höhlein et al., 2020). Recognizing this limitation, Shi et al. (2015) proposed a convolutional LSTM (ConvLSTM) architecture that combines the strengths of both LSTM networks and CNNs. In practical applications, the combination of LSTM networks and CNNs in the form of ConvLSTM models allowed for improvement of the accuracy of precipitation nowcasting (Shi et al., 2015), flood forecasting (Moishin et al., 2021), and river runoff forecasting (Ha et al., 2021; Zhu et al., 2023).

We use the Baltic Sea catchment as an example to illustrate our approach. Although the methodology we propose is universally applicable across various geographical regions, the Baltic Sea is a challenging region due to its unique hydrological characteristics, which are nearly decoupled from the open ocean (see Fig. 1) (Meier and Döscher, 2002). Freshwater enters the Baltic Sea through river runoff or positive net precipitation (precipitation minus evaporation) over the sea surface. The net precipitation accounts for 11 % and the river input for 89 % of the total freshwater input (Meier and Döscher, 2002). Consequently, the Baltic's sea surface salinity (SSS) is largely driven by freshwater supply from rivers.

Figure 1 Map of the Baltic Sea region. The red dots represent the locations of the major rivers that flow into the Baltic Sea, as represented in the E-HYPE hydrological model. The annotation BY15 marks the chosen validation station situated in the central Baltic Sea, which is used for validating the regional ocean model. The four remaining annotations (red) indicate the positions of the rivers that will be evaluated in detail.

SSS is an essential variable for the marine ecosystem in the Baltic Sea, as most species are adapted to either marine or freshwater conditions. Therefore, biases in the SSS's spatial and temporal variability significantly impact primary production and fish biomass (Kniebusch et al., 2019). Accurate modeling of the Baltic Sea relies heavily on the quality of the river input data used in the simulations. Analyzing nearly 100 years of observations, Winsor et al. (2001) found that variations in freshwater storage are closely correlated with accumulated changes in river runoff. From 1902 to 1998, the average freshwater inflow into the Baltic amounted to 16 115 m³ s⁻¹, with contributions from river runoff (14 085 m³ s⁻¹) and net precipitation over the Baltic Sea (2030 m³ s⁻¹) (Meier and Kauker, 2003). This freshwater inflow results in a residence time of about 35 years for freshwater in the Baltic Sea (Omstedt and Hansson, 2006; Winsor et al., 2001; Meier and Kauker, 2003).

In this work, we demonstrate that ConvLSTM networks are a reliable method for predicting multiple rivers simultaneously, using only atmospheric forcing as input data, which allows us to emulate a hydrological model. The main focus of this work is on presenting a ConvLSTM architecture capable of predicting daily river runoff for 97 rivers across the Baltic Sea catchment. To train the network, we use data from the E-HYPE model (Väli et al., 2019) as reference output data and data from the UERRA project (Uncertainties in Ensembles of Regional Reanalyses, http://www.uerra.eu/, last access: 25 March 2025) as atmospheric forcing (Sect. 3). The quality of the model is evaluated in Sect. 4. The obtained results are discussed further and evaluated in Sect. 5.

2 Implemented model architecture

2.1 The main idea

We assume that the runoff at a specific point in time t for all N_r considered rivers collected in the vector ${\mathit{R}}^t\in {\mathbb{R}}^{N_r}$ can be accurately approximated by a functional $\mathit{M}\left(\mathit{\right\{}{X}^{t_k}[x,y,\mathit{\tau }]\mathit{\left\}}\right)$ of $k=\mathrm{1},\mathrm{\dots },N_k$ atmospheric fields $X^{t_k}[x,y,\mathit{\tau }]$ which are known for the past $\mathit{\tau }=\mathrm{1},\mathrm{\dots },N_{\mathit{\tau }}$ time instances. This relationship is expressed as

$$\begin{array}{}\text{(1)}& {\mathit{R}}^t=\mathit{M}\left(\mathit{\right\{}{X}_k^t[x,y,\mathit{\tau }]\mathit{\left\}}\right).\end{array}$$

The atmospheric fields are evaluated over a spatial domain $x=\mathrm{1},\mathrm{\dots },N_x$ and $y=\mathrm{1},\mathrm{\dots },N_y$, which is sufficiently large to capture all significant local and non-local contributions of the atmospheric fields to the river runoff. Typically, such mapping is realized using a hydrological model that simulates all relevant physical processes, transforming variables like precipitation and evaporation into river runoff. This process relies heavily on domain knowledge to tune all parameters to reasonable values. As an alternative, combining a ConvLSTM model with a subsequent fully connected (FC) neural network can adequately represent this functional. This approach eliminates the need for detailed knowledge of the involved processes and their modeling. Instead, these features can be “learned” by the network in an automated manner; i.e., all free parameters are optimized such that the network's output reproduces the data of the hydrological model with given atmospheric input fields. Our proposed network architecture is visualized in Fig. 2 and is described in detail in the following sections. To provide an overview, we will discuss the main components of this architecture one by one.

Figure 2 Combined ConvLSTM and FC network architecture. The starting point is the continuous time series of input data {X^t} (upper white block). From this series, a contiguous sequence of N_τ elements (yellow block) is used to feed a chain of N_τ connected ConvLSTM cells (light-blue block) building the ConvLSTM network (grey block). The input sequence is mapped via weighting matrices 𝓜^g (green blocks) onto gate vectors g^t[τ]. The gate vectors are then used to update the cell state C^t[τ−1] and the hidden state H^t[τ−1] of the last ConvLSTM cell to the current values C^t[τ] and H^t[τ], respectively. The update is performed with the LSTM core equations collectively described by the mapping L; see Eq. (4). The weighting matrices 𝓝^g (green blocks) control how much of the last hidden state enters the updated state. The final output of the ConvLSTM H^t[τ] is then propagated to a FC network, which itself is a chain of three FC layers consisting of weighting matrices 𝓦 and connected via ReLU functions; see Sect. 2.4. The final result is the river runoff R^t for all rivers considered at the current time instance t (white block in the lower left corner). Note that all the bias vectors are omitted for the sake of clarity. See the text for more information.

2.2 The ConvLSTM network

2.2.1 The LSTM approach

Before turning directly to the ConvLSTM, the simpler architecture of the plain LSTM model is examined. This serves as a foundation for understanding the more complex ConvLSTM. The LSTM network, a specialized form of a recurrent neural network (RNN), is specifically designed to model temporal sequences ${\mathit{X}}^t\left[\mathrm{1}\right],\mathrm{\dots }{\mathit{X}}^t\left[\mathit{\tau }\right],\mathrm{\dots }{\mathit{X}}^t\left[N_{\mathit{\tau }}\right]$ of N_τ input quantities ${\mathit{X}}^t\left[\mathit{\tau }\right]=\left(X_k^t\right[\mathit{\tau }\left]\right)\in {\mathbb{R}}^{N_k}$. This sequence is taken from a dataset given in the form of a time series {X^t}, with the endpoint coinciding with the specific element in the time series connected to time t, i.e., X^t[N_τ]≡X^t; see Fig. 2. Here, N_k represents the number of input “channels”, which can correspond to different measurable quantities. The LSTM's unique design allows it to adeptly handle long-range dependencies, setting it apart from traditional RNNs in terms of accuracy (see Fig. 3).

Figure 3 Inner structure of a long short-term memory cell. See Fig. 2 and the text for information.

The critical component of the LSTM's innovation is its cell state, ${\mathit{C}}^t\left[\mathit{\tau }\right]=\left(C_h^t\right[\mathit{\tau }\left]\right)\in {\mathbb{R}}^{N_h}$, which stores state information also referred to as long-term memory. This state information complements the so-called hidden state vector ${\mathit{H}}^t\left[\mathit{\tau }\right]=\left(H_h^t\right[\mathit{\tau }\left]\right)\in {\mathbb{R}}^{N_h}$, which is also known from simpler neural network architectures. In the case of the LSTM, the hidden state vector plays the role of short-term memory. The cell state and the hidden state are vectors where each element is associated with one of the N_h hidden layers labeled by h. These internal, artificial degrees of freedom enable the high adaptability of neural networks. The two state vectors are determined through several self-parameterized gates, all in the same vector space as C^t[τ]; see Fig. 3 for a visualization.

In particular, the forget gate F^t[τ] defines the portion of the previous (long-term memory) cell state C^t[τ−1] that should be kept (see the dashed cyan box therein). The input gate I^t[τ] controls the contribution of the current input used to update the long-term memory C^t[τ] (magenta and yellow boxes). The output gate, O^t[τ], then determines how much this updated long-term memory contributes to the new (short-term memory) hidden state H^t[τ] (dashed black box).

For a fixed point τ in the sequence, the action of an LSTM cell, i.e., the connection between the input X^t[τ], the various gates, and the state vectors, is given mathematically as follows. First, the elements of the input sequence together with the hidden state are mapped onto auxiliary gate vectors, collectively denoted by ${\mathit{g}}^t\left[\mathit{\tau }\right]=\left(g_h^t\right[\mathit{\tau }\left]\right)\in {\mathbb{R}}^{N_h}$, via

$$\begin{array}{}\text{(2)}& g_h^t\left[\mathit{\tau }\right]={\mathcal{M}}_{hk}^g{X}_k^t\left[\mathit{\tau }\right]+{\mathcal{N}}_{h{h}^{\prime }}^g{H}_{h^{\prime }}^t[\mathit{\tau }-\mathrm{1}]+{\mathcal{B}}_h^g,\end{array}$$

where $g\in i,f,o,c$ stands for the input, forget, output, and cell state gates, respectively, and Einstein's summation convention is employed; i.e., indices that appear twice are summed over. The calligraphic symbols ${\mathcal{M}}_{hk}^g,{\mathcal{N}}_{h{h}^{\prime }}^g$, and ${\mathcal{B}}_h^g$ are the free parameters of the network that are optimized for the given problem during the training, which is at the heart of any machine learning approach. The matrix ${\mathcal{M}}^g=\left({\mathcal{M}}_{hk}^g\right)\in {\mathbb{R}}^{N_h\times N_k}$ can be interpreted as a Markovian-like contribution of the current input X^t[τ] to the gates, whereas ${\mathcal{N}}^g=\left({\mathcal{N}}_{h{h}^{\prime }}^g\right)\in {\mathbb{R}}^{N_h\times N_h}$ scales a non-Markovian part determined by the hidden state of the last sequence point τ−1. The vector ${\mathcal{B}}^g=\left({\mathcal{B}}_h^g\right)\in {\mathbb{R}}^{N_h}$ is a learnable bias. It should be stressed that these parameters do not depend on t or τ and are thus optimized once for the complete dataset {X^t}.

Note that this mapping is sometimes extended by a contribution to $g_h^t\left[\mathit{\tau }\right]$ from the past cell state C^t[τ−1]. Nevertheless, this “peeping” mechanism is not considered further in this work. For the sake of brevity, we can write the mapping more compactly in matrix vector form as

$$\begin{array}{}\text{(3)}& {\mathit{g}}^t\left[\mathit{\tau }\right]={\mathcal{M}}^g{\mathit{X}}^t\left[\mathit{\tau }\right]+{\mathcal{N}}^g{\mathit{H}}^t[\mathit{\tau }-\mathrm{1}]+{\mathcal{B}}^g.\end{array}$$

Second, the actual gate vectors are computed by the core equations of the LSTM as proposed by Hochreiter and Schmidhuber (1997):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{I}}^t\left[\mathit{\tau }\right]& =\mathit{\sigma }\left({\mathit{i}}^t\right[\mathit{\tau }\left]\right),\\ {\mathit{F}}^t\left[\mathit{\tau }\right]& =\mathit{\sigma }\left({\mathit{f}}^t\right[\mathit{\tau }\left]\right),\\ {\mathit{O}}^t\left[\mathit{\tau }\right]& =\mathit{\sigma }\left({\mathit{o}}^t\right[\mathit{\tau }\left]\right),\\ {\mathit{C}}^t\left[\mathit{\tau }\right]& ={\mathit{F}}^t\left[\mathit{\tau }\right]\circ {\mathit{C}}^t[\mathit{\tau }-\mathrm{1}]+{\mathit{I}}^t\left[\mathit{\tau }\right]\circ \tanh \left({\mathit{c}}^t\right[\mathit{\tau }\left]\right),\\ {\mathit{H}}^t\left[\mathit{\tau }\right]& ={\mathit{O}}^t\left[\mathit{\tau }\right]\circ \tanh \left({\mathit{C}}^t\right[\mathit{\tau }\left]\right),\end{array}\end{array}$$

where σ denotes the logistic sigmoid function, tanh is the hyperbolic tangent, and ∘ stands for the Hadamard product (all applied in an element-wise fashion to the vectors). In the last two equations, the roles of the input, forget, and output gates as described above become apparent.

The third step in a single-layer LSTM (as employed for the work presented here) is to provide the outputs of the current LSTM cell, i.e., H^t[τ] and C^t[τ], to the subsequent LSTM cell that processes the next element X^t[τ+1] of the input sequence.

The full action of the LSTM network up to the end of the sequence can be written as a nested function call:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \left({\mathit{H}}^t\left[N_{\mathit{\tau }}\right],{\mathit{C}}^t\left[N_{\mathit{\tau }}\right]\right)\\ & =L\left({\mathit{X}}^t\right[N_{\mathit{\tau }}],L({\mathit{X}}^t[N_{\mathit{\tau }}-\mathrm{1}],\mathrm{\dots }\\ & L\left({\mathit{X}}^t\left[\mathrm{1}\right],\left({\mathit{H}}^t\right[\mathrm{0}],{\mathit{C}}^t[\mathrm{0}\left]\right)\right)\mathrm{\dots }\left)\right),\end{array}\end{array}$$

where $L\left({\mathit{X}}^t\left[\mathit{\tau }\right],\left({\mathit{H}}^t[\mathit{\tau }-\mathrm{1}],{\mathit{C}}^t[\mathit{\tau }-\mathrm{1}]\right)\right)$ represents Eqs. (3) and (4). For the present work, the initial conditions are chosen as ${\mathit{H}}^t\left[\mathrm{0}\right]={\mathit{C}}^t\left[\mathrm{0}\right]=\mathrm{0}$, which means that there is no memory longer than N_τ time steps.

The final outputs of the ConvLSTM chain, H^t[N_τ] and C^t[N_τ], encode information on the entire input sequence, ending at time t. This information must be decoded via an appropriate subsequent network, as described in Sect. 2.4.

2.3 Combining the LSTM with spatial convolution

Although the plain LSTM performs well in handling temporal sequences of point-like quantities, it is not designed to recognize spatial features in sequences of two-dimensional maps as atmosphere–ocean interface fields. To address this limitation, we employ a ConvLSTM architecture as described below.

In this type of network, the elements of the input sequence are given as spatially varying fields ${\mathit{X}}^t\left[\mathit{\tau }\right]=\left(X_k^t\right[x,y,\mathit{\tau }\left]\right)\in {\mathbb{R}}^{N_k\times (N_x\times N_y)}$, where $x\in [\mathrm{1},N_x]$ and $y\in [\mathrm{1},N_y]$ run over the horizontal and vertical dimensions of the map. To enable the “learning” of spatial patterns, the free parameters of the network are replaced with two-dimensional convolution kernels ${\mathcal{M}}^g=\left({\mathcal{M}}_{hk}^g\right[\mathit{\xi },\mathit{\eta }\left]\right)\in {\mathbb{R}}^{(N_h\times N_k)\times (N_{\mathit{\xi }}\times N_{\mathit{\eta }})}$ and ${\mathcal{N}}^g=\left({\mathcal{N}}_{h{h}^{\prime }}^g\right[\mathit{\xi },\mathit{\eta }\left]\right)\in {\mathbb{R}}^{(N_h\times N_h)\times (N_{\mathit{\xi }}\times N_{\mathit{\eta }})}$. The width and height of the kernels are given by N_ξ and N_η, respectively, and $\mathit{\xi }\in [-(N_{\mathit{\xi }}-\mathrm{1})/\mathrm{2},(N_{\mathit{\xi }}-\mathrm{1})/\mathrm{2}],\mathit{\eta }\in [-(N_{\mathit{\eta }}-\mathrm{1})/\mathrm{2},(N_{\mathit{\eta }}-\mathrm{1})/\mathrm{2}]$. Without loss of generality, we assume odd numbers for the kernel sizes.

A convolution with these kernels then gives the mapping from the input quantities to the gates:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{g}_h^t[x,y,\mathit{\tau }]& ={\mathcal{M}}_{hk}^g[\mathit{\xi },\mathit{\eta }]X_k^t[x-\mathit{\xi },y-\mathit{\eta },\mathit{\tau }]\\ & +{\mathcal{N}}_{h{h}^{\prime }}^g[\mathit{\xi },\mathit{\eta }]H_{h^{\prime }}^t[x-\mathit{\xi },y-\mathit{\eta },\mathit{\tau }-\mathrm{1}]+{\mathcal{B}}_h^g,\end{array}\end{array}$$

again with Einstein's convention imposed.

It immediately becomes apparent that, in the case of the ConvLSTM, the gate and state vectors must become vector fields ($\in {\mathbb{R}}^{N_h\times (N_x\times N_y)}$) as well. We can write this mapping in the same way as Eq. (3) but by replacing the standard matrix vector multiplication with a convolution (denoted with ∗), i.e.,

$$\begin{array}{}\text{(7)}& {\mathit{g}}^t\left[\mathit{\tau }\right]={\mathcal{M}}^g\ast {\mathit{X}}^t\left[\mathit{\tau }\right]+{\mathcal{N}}^g\ast {\mathit{H}}^t[\mathit{\tau }-\mathrm{1}]+{\mathcal{B}}^g.\end{array}$$

The subsequent processing of g^t[τ] remains symbolically the same as presented in Eq. (4) but with all appearing quantities now meaning vector fields instead of simple vectors.

In summary, the ConvLSTM is designed to process tasks that require a combined understanding of spatial patterns and temporal sequences in data. It merges the image-processing capabilities of CNNs with the time series modeling of LSTM networks.

2.4 Fully connected layer

As stated in Sect. 2.2.1, the final outputs H^t[N_τ] and C^t[N_τ] of the ConvLSTM encode information on the full input sequence. To contract this information to obtain the runoff vector R^t representing the N_r rivers, we propose subjecting the final short-term memory (i.e., the hidden state H^t[N_τ]) to an additional FC network.

In particular, the dimensionality of the vector field H^t[N_τ] is reduced sequentially by three nested FC layers, each connected to the other by the rectified linear unit (ReLU); see Fig. 2. Integrating out artificial degrees of freedom in a step-wise fashion has turned out to be beneficial.

The runoff of the rth river is then obtained via (using Einstein's convention)

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{R}_r^t& ={\mathcal{W}}_{rb}^{\mathrm{3}}\mathrm{ReLU}\left({\mathcal{W}}_{ba}^{\mathrm{2}}\mathrm{ReLU}\right({\mathcal{W}}_{ah}^{\mathrm{1}}[x,y]H_h^t[x,y,N_{\mathit{\tau }}]\\ & +{\mathcal{B}}_a^{\mathrm{1}})+{\mathcal{B}}_b^{\mathrm{2}})+{\mathcal{B}}_r^{\mathrm{3}},\end{array}\end{array}$$

where $a=\mathrm{1},\mathrm{\dots },N_a$, $b=\mathrm{1},\mathrm{\dots },N_b$, and the hyperparameters N_a and N_b are chosen such that $N_h\cdot N_x\cdot N_y>N_a>N_b>N_r$ in order to achieve the aforementioned step-by-step reduction in dimensionality. The weights 𝒲 and biases ℬ stand for parameters that are optimized during the training of the network.

In matrix vector notation, this can be compressed to

$$\begin{array}{}\text{(9)}& {\mathit{R}}^t={\mathcal{W}}^{\mathrm{3}}\mathrm{ReLU}\left({\mathcal{W}}^{\mathrm{2}}\mathrm{ReLU}\left({\mathcal{W}}^{\mathrm{1}}{\mathit{H}}^t\left[N_{\mathit{\tau }}\right]+{\mathcal{B}}^{\mathrm{1}}\right)+{\mathcal{B}}^{\mathrm{2}}\right)+{\mathcal{B}}^{\mathrm{3}}.\end{array}$$

Combining Eq. (9) with Eq. (5) finally provides an explicit formula for the initial assumption of modeling the runoff for time t as a function of a sequence of atmospheric fields; i.e.,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{R}}^t& =\mathit{M}\left(\mathit{\right\{}{X}_k^t[x,y,\mathit{\tau }]\mathit{\left\}}\right)\\ & ={\mathcal{W}}^{\mathrm{3}}\mathrm{ReLU}\left({\mathcal{W}}^{\mathrm{2}}\mathrm{ReLU}\right({\mathcal{W}}^{\mathrm{1}}{\mathit{L}}_{\mathrm{H}}\left({\mathit{X}}^t\right[N_{\mathit{\tau }}],\\ & L\left({\mathit{X}}^t[N_{\mathit{\tau }}-\mathrm{1}],\mathrm{\dots }L\left({\mathit{X}}^t\left[\mathrm{1}\right],(\mathrm{0},\mathrm{0})\right)\mathrm{\dots }\right))+{\mathcal{B}}^{\mathrm{1}})\\ & +{\mathcal{B}}^{\mathrm{2}})+{\mathcal{B}}^{\mathrm{3}},\end{array}\end{array}$$

where L_H means that only the hidden state vector of the final ConvLSTM call is forwarded to the FC layer.

3 Technical details

3.1 Runoff data used for training

The non-stationary daily runoff data covering the period 1979 to 2011 are based on an E-HYPE hindcast simulation that was forced by a regional downscaling of ERA-Interim (Dee et al., 2011) with RCA3 (Samuelsson et al., 2011) and implemented in NEMO-Nordic (Hordoir et al., 2019) as a mass flux. The BMIP project (Gröger et al., 2022) played a crucial role in addressing the lack of consistent river discharge data for the entire study period (1961–2018). On this point, no comparable long-term dataset with daily resolution was available. In other studies multiple datasets have been merged but offer only monthly resolution (see, e.g., Fig. 3, Meier et al., 2019). Hence, a new homogeneous runoff dataset was created. The 1961–1978 runoff data are based on Bergström and Carlsson (1994), with values interpolated from monthly to daily scales. The 2012–2018 data are derived from an E-HYPE forecast product. To ensure consistency for the analysis, the periods before (1961 to 1978) and after (2012 to 2018) have been neglected. Notably, the Neva River is an exception, as its discharge data originate from observational records (1961–2016) provided by the Russian State Hydrological Institute rather than E-HYPE hindcasts.

It should be noted that, for this study, we used an intermediate dataset of river runoff developed during BMIP that was employed to run the ocean model. In this dataset, some rivers had not yet been merged, resulting in discrepancies between the totals of freshwater input locations of 97 in this study and 91 rivers in the final version of Väli et al. (2019). The quality of the runoff was extensively evaluated. The dataset was found to closely align with historical observations for various rivers and with the Bergström and Carlsson (1994) dataset, showing a difference of under 1 % for the total Baltic Sea runoff (Väli et al., 2019). For more information, see Gröger et al. (2022) and Väli et al. (2019).

3.2 Atmospheric forcing

The UERRA-HARMONIE regional reanalysis dataset was developed as part of the FP7 UERRA project (http://www.uerra.eu/). The UERRA-HARMONIE system represents a comprehensive, high-resolution reanalysis covering a wide range of essential climate variables. This dataset includes data on air temperature, pressure, humidity, wind speed and direction, cloud cover, precipitation, albedo, surface heat fluxes, and radiation fluxes from January 1961 to July 2019. With a horizontal resolution of 11 km and analyses conducted at 00:00, 06:00, 12:00, and 18:00 UTC, it also provides hourly resolution forecast model data. UERRA-HARMONIE is accessible through the Copernicus Climate Data Store (CDS, https://cds.climate.copernicus.eu/#!/home, last access: 25 March 2025) initially produced during the UERRA project and later changed to the Copernicus Climate Change Service (C3S, https://climate.copernicus.eu/copernicus-regional-reanalysis-europe, last access: 25 March 2025). For the training of the neural network, the hourly data were remapped to daily values.

Lastly, it should be noted that UERRA is not the atmospheric dataset that was used to drive the original E-HYPE model.

3.3 Ocean model

We use a coupled three-dimensional ocean model, called the Modular Ocean Model (MOM) (Griffies, 2012), to simulate the Baltic Sea. It has a horizontal resolution of 3 nautical miles, roughly corresponding to 5.556 km, and 152 vertical z^* levels with a first layer thickness of 0.5 m and a total depth of 264 m. This model uses a finite-difference method to solve the full set of primitive equations to calculate the motion of water and the transport of heat and salt. The K-profile parameterization (KPP) was used as a turbulence closure scheme. The model's western boundary opens into the Skagerrak and connects the Baltic Sea to the North Sea. A more detailed description of the setup can be found in Gröger et al. (2022).

3.4 Neural network hyperparameters

Our architecture is implemented as a sequential model, which allows for testing of multiple ConvLSTM layers – a concatenation of multiple ConvLSTM cells. The following ConvLSTM output is then mapped by three fully connected linear layers, where the final layers map the output to 97 rivers. We used a custom loss function similar to a mean squared error loss that penalizes outliers more severely. For the training, we use the AdamW optimizer, which is an improved version of Adam that decouples weight decay (regularization) from the gradient update, likely leading to better generalization. The optimizer is configured with a learning rate and a weight decay to prevent overfitting. To adapt the learning rate during training, we employ the ReduceLROnPlateau scheduler. This scheduler monitors the validation mean squared error (MSE) and reduces the learning rate by a factor of 10 if no improvement is observed over 10 epochs. This dynamic adjustment helped the model converge more efficiently and avoid getting stuck in local minima.

The best set of hyperparameters has been defined by iterating over a pre-defined selection of possible parameters. The set of hyperparameters that has been chosen for the present study is given in Table 1. The model's performance can be described as relatively robust when changing the set of hyperparameters (see Fig. A3). Interestingly, smaller input sizes of 10 d also perform really well. However, we still decided to use longer timescales, as we assume that larger input sizes increase the stability of the model needed for long-term climate simulations.

Table 1 Final set of parameters for the ConvLSTM model.

River runoff is influenced not only by the current day's atmospheric conditions, but also by the cumulative and lagged effects of prior days' weather patterns. The choice of atmospheric fields was based on the assumption that the runoff should be mainly given by the net precipitation (precipitation minus evaporation). The evaporation flux is often calculated as a function of wind speed, the air's humidity and density, and the involved turbulent exchange coefficients (Karsten et al., 2024), where the air temperature influences the latter two. Hence, as input, the model receives N_τ=30 d of the atmospheric surface fields temperature, precipitation, specific humidity, and wind speed, with a daily resolution to predict daily river runoff R (see Fig. 4). This window size allows the model to “remember” key atmospheric conditions leading up to a given day, enabling it to accurately predict runoff.

Figure 4 Schematic structure of the ConvLSTM implementation for river runoff forecasting.

4 Results

4.1 ConvLSTM model evaluation

The model was trained and evaluated with daily data from 1979 to 2011, as this period represents the only period of E-HYPE without further bias correction applied to the runoff to match the observations. The complete dataset was divided into randomly chosen splits of 80 % training data, 10 % validation data to evaluate the model's performance during training, and 10 % test data which are finally used to assess the model's performance after training. The model was trained for 400 epochs, and the model weights with the lowest mean squared validation error have been stored. The model's accuracy for the combined daily predicted runoff from all 97 rivers flowing into the Baltic Sea is displayed in Fig. 5. For this evaluation, the model's output is compared to the test data, which the model has not seen during the training phase.

Panel (a) illustrates the relative prediction error in relation to the E-HYPE data, indicating that, on daily timescales, the model can predict river runoff with an accuracy of ±5 %. The overall correlation is 0.997, with the resulting error metrics yielding a root mean square error (RMSE) of 323.99 m³ s⁻¹ and a mean absolute error (MAE) of 249.51 m³ s⁻¹. While the model's performance is already satisfactory, the discrepancies between the actual values and the predictions can partly be attributed to the use of a different atmospheric dataset than the one originally used to drive the E-HYPE model. However, by applying a rolling mean with a 5 d window, the prediction error is reduced to less than 1 %, which is acceptable for climate modeling purposes. Panel (b) displays the distribution of residuals as a density plot. The figure shows that the distribution of residuals follows a Gaussian shape. The bell-shaped curve is approximately centered around zero, indicating that the model does not exhibit a systematic bias and meaning that it does not consistently overestimate or underestimate the river runoff values. Most residuals lie within a narrow range around zero, suggesting that large prediction errors are relatively rare.

For individual rivers, the distribution of the residuals still follows a Gaussian shape. However, on daily timescales the errors are larger, reaching ±30 % during individual peaks (Fig. A4 in the Appendix).

Figure 5 Model accuracy for predicted river runoff. (a) Relative prediction error over time in relation to the E-HYPE data. (b) Density plot of the residuals.

In the following, the model's performance in reproducing the total river runoff and the discharges of four individual rivers into the Baltic Sea is addressed. Using the test dataset, Fig. 6 shows the predicted river runoff and E-HYPE data.

Panel (a) illustrates the total river runoff into the Baltic Sea, with both the predicted runoff (ConvLSTM) and the river runoff of the hydrological model E-HYPE smoothed using a 5 d rolling mean. The predicted river runoff closely matches the original data, demonstrating the model's accuracy in predicting the overall river runoff into the Baltic Sea. Panel (b) focuses on the Neva River (see Fig. 1), one of the largest rivers flowing into the Baltic Sea. The residual plot illustrates the prediction errors relative to the E-HYPE river runoff data over time. The ConvLSTM model predicts the Neva River runoff within a ± 2.5 % range.

Panels (c)–(e) show the prediction residuals for the Oder River, the Ume River, and the Neman River (see Fig. 1). Compared to the Neva River, the prediction errors are larger, which may be attributed to the training dataset, as the runoff for the Neva River is based on measurements, whereas that of the other rivers is based solely on E-HYPE. Still, all of the results lie within the error margin of E-HYPE itself compared to the observations (Fig. A1), with the average error on daily timescales for individual rivers mostly under 10 %, showcasing the model's ability to forecast runoff for this river accurately.

The residuals were calculated as the relative difference between the predicted runoff and runoff data used for training and were finally normalized by the runoff data used for training.

An overview of the individual error metrics is given in Table A1.

Figure 6 Model performance in predicting river runoff: (a) total river runoff to the Baltic Sea with a 5 d rolling mean for both the predicted and original data of the E-HYPE hydrological model. (b–e) Residuals of runoff prediction for individual rivers showing the prediction error over time. The residuals were calculated as the relative difference between the predicted and observed values, normalized by the observed values.

4.2 Application of the ConvLSTM in combination with an ocean model

Lastly, we evaluate the performance of the ConvLSTM by incorporating the predicted river runoff as hydrological forcing into MOM5. This provides a robust validation of the runoff model against more complex real-world conditions and ensures that the predictions accurately reflect the impact of the river discharge on the ocean dynamics. This, in turn, confirms that the ConvLSTM captures the temporal and spatial variability of river runoff and that the residuals shown in Fig. 5 are indeed insignificant when it comes to realistic applications.

Figure 7 shows the salinity comparison between E-HYPE and the predicted river runoff at BY15 – a central station in the Baltic Sea east of Gotland. The ocean model simulation using the river runoff predicted by the ConvLSTM closely mirrors the control simulation that is forced with the E-HYPE runoff. Panel (a) shows the surface salinity, representing the high-frequency variations in salinity heavily affected by river runoff. The predicted salinity using ConvLSTM river runoff matches the control simulation well, capturing the short-term fluctuations effectively. Panel (b) shows the bottom salinity, representing low-frequency variations in the Baltic Sea, which is also reproduced well with the ConvLSTM predictions. It should be noted that the discrepancies between the simulated salinity and the observed values at BY15 are not directly linked to the performance of the ConvLSTM river runoff model. Instead, they are attributed to MOM5's representation of physical processes, particularly the treatment of mixing, advection, and stratification in the Baltic Sea. Several factors may contribute to these discrepancies. The Baltic Sea is known for its strong vertical stratification due to the input of freshwater from rivers. MOM5 uses the KPP scheme for turbulence, which may not perfectly resolve small-scale mixing processes and vertical salinity gradients. This can result in an overestimation of salinity variability at the surface. Moreover, while MOM5 captures the large-scale dynamics of the Baltic Sea, the lateral transport of saltwater from the Skagerrak into the central Baltic Sea may not be represented perfectly. This can introduce variability in surface and bottom salinity that is not observed in reality. However, all in all, the long-term trends and larger salinity changes are captured accurately, indicating the model's robustness in predicting high-frequency and low-frequency variations.

Figure 7 Salinity comparison at BY15 in the Baltic Sea: (a) surface salinity showing high-frequency variations and (b) bottom salinity showing low-frequency trends. The ConvLSTM predictions follow the E-HYPE data closely, demonstrating the model's accuracy in reproducing salinity levels affected by river runoff.

The final evaluation of the ConvLSTM model focuses on the spatial accuracy of river runoff predictions as visualized in Fig. 8. Panel (a) exhibits the vertically averaged salinity from 1981 to 2011 in the reference simulation. It highlights the Baltic Sea's strong horizontal gradients and complex topographic features, as evidenced by salinity variations in deeper waters captured by the vertical integration. In panel (b), these reference results are compared to the ConvLSTM simulation by showing the percentage difference in vertically averaged salinity. Overall, the differences remain below 1 %, except in the Gulf of Riga (22–24° E, 56.5–58.5° N for orientation), where the Daugava River dominates the runoff. The difference is approximately 1 %.

Figure 8 Spatial accuracy of river runoff predictions using the ConvLSTM model. (a) Vertically averaged salinity for the period 1981 to 2011 in the reference simulation, highlighting the sharp horizontal gradients and complex topographic features in the Baltic Sea. (b) Percentage difference in vertically averaged salinity between the reference simulation and the ConvLSTM simulation.

5 Conclusions

With the increasing demand from decision-makers for regional climate projections to quantify regional climate change impacts, the availability of precise hydrological forecasting becomes invaluable. In this work, we describe the implementation of a ConvLSTM network for predicting river runoff in a regional climate model, highlighting its potential to enhance river runoff forecasting across different coastal seas. Our model not only reproduces the total river runoff entering the Baltic Sea but also performs exceptionally well for individual rivers.

All of the results lie within the error margin of the hydrological model itself when compared to the observations, with the average error on daily timescales for individual rivers mostly under 10 %. Hence, our results confirm the excellent performance of LSTM networks in predicting river runoff (Humphrey et al., 2016; Huang et al., 2014; Ashrafi et al., 2017; Liu et al., 2020; Fang and Shao, 2022; Kratzert et al., 2018). In addition, our results align well with the observed performance of ConvLSTM networks in similar applications for predicting individual rivers (Ha et al., 2021), basin-wide runoff (Zhu et al., 2023), and precipitation nowcasting (Shi et al., 2015). In our study, we further extend the use of ConvLSTM networks by predicting multiple (n=97) rivers at once while maintaining high accuracy for the entire Baltic Sea as well as individual rivers. Moreover, the predicted river runoff proved robust when using the river runoff in a comprehensive ocean model setup of the Baltic Sea. Extending the simulation beyond the known period also provided robust results (Fig. A2).

The transition from traditional hydrological models to machine learning approaches, such as the implementation of the ConvLSTM model, offers significant advantages as the model can be integrated seamlessly into regional climate models, allowing for real-time computation of river runoff while making climate projections. While the initial training of the model requires substantial computational resources, this remains significantly less intensive than running comprehensive hydrological models. Furthermore, once trained, the ConvLSTM model is computationally efficient during inference, ensuring enhanced forecasting capabilities without significantly increasing computational demands. The achieved speedup (depending on the complexity of the hydrological model) is within the range of 30 to 90 times faster.

Nevertheless, the quality of the ConvLSTM model still depends on the performance of the hydrological model, which provides a comprehensive, homogeneous dataset that is essential for effective training. While different in their architectures, both the hydrological model and the ConvLSTM model provide a mapping from atmospheric variables to river runoff. This training approach contrasts with using measurement data for training, which is significantly more challenging due to the data sparsity over larger regions and potentially varying measurement techniques. Thus, rather than rendering traditional methods obsolete, the integration of machine learning models builds upon and enhances the foundational data provided by hydrological model methods.

The robust performance of the ConvLSTM model in simulating river runoff and its possible effective integration into coupled regional climate models, as in the IOW ESM (Karsten et al., 2024), pave the way for a multitude of new storyline simulations. Importantly, this can now be achieved without any expert domain knowledge on hydrological modeling. Hence, we can now explore various “what if” scenarios more reliably, under the assumption that the model weights attained during the training are robust. Such scenario testing is crucial for crafting effective water resource management strategies and adapting to a changing climate, and it therefore represents a significant step forward in our ability to understand and predict the complex dynamics of river systems and their impact on regional climate systems.

Appendix A

Figure A1 Residual for the hydrological model E-HYPE as well as for the prediction by the ConvLSTM model. The residuals were calculated as the relative differences between the predicted and observed values, normalized by the observed values.

Figure A2 Extended regional ocean model simulation for the period 1961–2009. Station BY15 is validated.

Figure A3 Model accuracy for predicted river runoff for 10 different hyperparameters. (a) Relative prediction error over time in relation to the E-HYPE data. (b) Density plot of the residuals.

Table A1 Model performance in predicting river runoff, together with the error metrics.

Figure A4 Model accuracy for predicted river runoff for individual rivers. (a) Relative prediction error over time in relation to the E-HYPE data. (b) Density plots of the residuals.