# Advances in Land Surface Model-based Forecasting: A comparative study of
LSTM, Gradient Boosting, and Feedforward Neural Network Models as prognostic
state emulators

Marieke Wesselkamp
&Matthew Chantry
&Ewan Pinnington 22footnotemark: 2
&Margarita Choulga 22footnotemark: 2
&Souhail Boussetta 22footnotemark: 2
&Maria Kalweit
&Joschka Boedecker 33footnotemark: 3 444BrainLinks-BrainTools, Freiburg,
Germany
&Carsten F. Dormann 11footnotemark: 1
&Florian Pappenberger 22footnotemark: 2
&Gianpaolo Balsamo22footnotemark: 2 555World Meteorological Organization,
Geneva, Switzerland
Department of Biometry, University of Freiburg, GermanyEuropean Centre for
Medium-Range Weather Forecasts, Reading, United KingdomDepartment of Computer
Science, University of Freiburg, Germany

###### Abstract

Most useful weather prediction for the public is near the surface. The
processes that are most relevant for near-surface weather prediction are also
those that are most interactive and exhibit positive feedback or have key role
in energy partitioning. Land surface models (LSMs) consider these processes
together with surface heterogeneity and forecast water, carbon and energy
fluxes, and coupled with an atmospheric model provide boundary and initial
conditions. This numerical parametrization of atmospheric boundaries being
computationally expensive, statistical surrogate models are increasingly used
to accelerated progress in experimental research. We evaluated the efficiency
of three surrogate models in speeding up experimental research by simulating
land surface processes, which are integral to forecasting water, carbon, and
energy fluxes in coupled atmospheric models. Specifically, we compared the
performance of a Long-Short Term Memory (LSTM) encoder-decoder network,
extreme gradient boosting, and a feed-forward neural network within a physics-
informed multi-objective framework. This framework emulates key states of the
ECMWF’s Integrated Forecasting System (IFS) land surface scheme, ECLand,
across continental and global scales. Our findings indicate that while all
models on average demonstrate high accuracy over the forecast period, the LSTM
network excels in continental long-range predictions when carefully tuned, the
XGB scores consistently high across tasks and the MLP provides an excellent
implementation-time-accuracy trade-off. The runtime reduction achieved by the
emulators in comparison to the full numerical models are significant, offering
a faster, yet reliable alternative for conducting numerical experiments on
land surfaces.

##  1 Introduction

While forecasting of climate and weather system processes has long been a task
for numerical models, the recent development in deep learning has introduced
competitive machine-learning (ML) systems for numerical weather prediction
(NWP) [1, 2, 3]. Land surface models (LSMs), even though being an integral
part of numerical weather prediction, have not yet caught the attention of the
ML-community. LSMs forecast water, carbon and energy fluxes, and in coupling
with an atmospheric model, provide the lower boundary and initial conditions
[4, 5]. The parametrization of land surface states thus does not only affect
predictability of earth and climate systems on sub-seasonal scales [6], but
also the short- and medium-range skill of NWP forecasts [5]. Beyond the online
integration with NWPs, offline versions of LSMs provide research tools for
experiments on the land surface [4], the diversity of which are however
limited by the required substantial computational resources and often moderate
runtime efficiencies [7].

Emulators constitute statistical surrogates for numerical simulation models
that, by approximating the latter, aim at increasing computational efficiency
[8]. While for construction emulators can themselves require substantial
computational resources, their subsequent evaluation usually runs orders of
magnitude faster than the original numerical model [9]. For this reason,
emulators have found application for example in modular parametrization of
online weather forecasting systems [10], in replacing the MCMC-sampling
procedure in Bayesian calibration of ecosystem models [9], or in generating
ensembles of atmospheric states for forecast uncertainty quantification [11].
Beyond their computational efficiency, surrogate models with high parametric
flexibility have the potential to correct for process mis-specification and
improve predictions towards a physical model [12].

Modelling approaches used for emulation range from low parametrized, auto-
regressive linear models to highly non-linear and flexible neural networks
[13, 14, 15]. In the global land surface system M-MESMER, a set of simple AR1
regression models is used to initialize the numerical LSM, resulting in a
modularized emulator [13]. Numerical forecasts of gross primary productivity
and hydrological targets were successfully approximated by Gaussian processes
[8, 14], the advantage of which is their direct quantification of prediction
uncertainty. When it comes to highly diverse or structured data, neural
networks have shown to deliver accurate approximations for variables from
gravity wave drags to urban surface temperature [10, 15]. In most fields of
machine learning, specific types of neural networks are now the best approach
to representing fit and prediction. One exception is so-called tabular data,
i.e. data without spatial or temporal interdependencies (as opposed to vision
and sound), where extreme gradient boosting is still the go-to approach [16,
17].

ECLand is the land surface scheme that provides boundary and initial
conditions for the Integrated Forecasting System (IFS) of the European Centre
for Medium-range Weather Forecasts (ECMWF) [4]. Driven by meteorological
forcing and spatial climate fields, it has a strong influence on the NWP [5]
and also constitutes a standalone framework for offline forecasting of land
surface processes, the advantage of which for the online framework is the
temporal consistency of prognostic state variables [6]. The modular
construction of ECLand offers potential for element-wise improvement of
process representation and thus a stepwise development towards increased
computational efficiency. Within the IFS, ECLand also forms the basis of the
land surface data assimilation system, updating the land surface state with
synoptic data and satellite observations of soil moisture and snow. Emulators
of physical systems have been shown to be beneficial in data assimilation
routines, allowing for a quick and low maintenance estimation of the tangent
linear model [18, 5]. Together with the potential to run large ensembles of
land surface states at a much-reduced cost, this would be a potential
application of the surrogate models introduced here. Long-short term memory
networks (LSTMs) have gained popularity in hydrological forecasting as
rainfall-runoff models, for predicting stream flow temperature and also soil
moisture [e.g. [19, 20, 21, 22]]. Research on the interpretability of LSTMs
has found correlations between the model cell states and spatially or
thematically similar hydrological units [20], suggesting the specific
usefulness of LSTM for representing variables with dynamic storages and
reservoirs [23]. As emulators, LSTMs have been shown useful for sea surface
level projection in a variational manner with Monte Carlo dropout [24]. While
most of these studies trained their models on observations or reanalysis data,
our emulator learns the representation from ECLand simulations directly. To
our knowledge, a comparison of models without memory mechanisms to an LSTM-
based neural network for global land surface emulation has not been conducted
before.

We emulate seven prognostic state variables of ECLand, which represent core
land surface processes: soil water volume and soil temperature, each at three
depth layers, and snow cover fraction at the surface layer. These three state
variables represent the core of the current configuration of ECLand We
specifically focus on the utility of memory mechanisms, highlighting the
development of a single LSTM-based encoder-decoder model compared to an
extreme gradient boosting approach (XGB) and a multilayer perceptron (MLP),
which all perform the same tasks. The LSTM architecture builds on an encoder-
decoder network design introduced for flood forecasting [25]. To compare
forecast skill systematically, the three emulators were compared in short- and
long-range forecasting against climatology [26]. In this work, evaluation is
done on ECLand simulations only, i.e. on purely synthetic data, while future
work will encompass transfer learning and validation on observations.

##  2 Methods

###  2.1 The Land Surface Model: ECLand

ECLand is a tiled ECMWF Scheme for Surface Exchanges over Land that represents
surface heterogeneity and incorporates land surface hydrology (ECLand) [27,
28]. ECLand computes surface turbulent fluxes (of heat, moisture and momentum)
and skin temperature over different tiles (vegetation, bare soil, snow,
interception and water) and then calculates an area-weighted average for the
grid-box to couple with the atmosphere [4]. For the overall accuracy of the
model, accurate parameterizations are essential [29] as e.g. the land surface
parameterization determines the sensible and latent heat fluxes, and provide
the lower boundary conditions for enthalpy and moisture equations in the
atmosphere [30]. We emulate three prognostic state variables of ECLand, that
represent core land surface processes: soil water volume and soil temperature
at each three depth layers (each at 0 – 7 cm, 7 – 21 cm and 21 – 72 cm) and
snow cover fraction, aggregated at the surface layer, so below are some more
details on these parametrisations.

###  2.2 Data sources

As training data base, global simulation and reanalysis time series from 2010
to 2022 were compiled to zarr format at an aggregated 6-hourly temporal
resolution. Simulations and climate fields were generated from ECMWFs
development cycle CY49R2, ECland forced by ERA-5 meteorological reanalysis
data [31]. There are three main sources of data used for creation of the data
base: The first is a selection of surface physiographic fields from ERA5 [31]
and their updated versions [32, 4, 6] used as static model input features (X).
The second is a selection of atmospheric and surface model fields from ERA5,
used as static and dynamic model input features (Y). The third is ECLand
simulation results, constituting the model’s dynamic prognostic state
variables (z) and hence model input and target features. A total of 41 static,
seasonal and dynamical features were used to create the emulators, see Table 1
for an overview of input variables and details on the surface physiographic
and atmospheric fields below.

####  2.2.1 Surface physiographic fields

Surface physiographic fields have gridded information of the Earth’s surface
properties (e.g. land use, vegetation type, and distribution) and represent
surface heterogeneity in the ECLand of the IFS [29]. They are used to compute
surface turbulent fluxes (of heat, moisture, and momentum) and skin
temperature over different surfaces (vegetation, bare soil, snow,
interception, and water) and then to calculate an area-weighted average for
the grid box to couple with the atmosphere. To trigger all different
parametrization schemes, the ECMWF model uses a set of physiographic fields
that do not depend on initial condition of each forecast run or the forecast
step. Most fields are constant; surface albedo is specified for 12 months to
describe the seasonal cycle. Depending on the origin, initial data come at
different resolutions and different projections and are then first converted
to a regular latitude–longitude grid (EPSG:4326) at approx. 1km at Equator
resolution and secondly to a required grid and resolution. Surface
physiographic fields used in this work consist of orographic, land, water,
vegetation, soil, albedo fields, see Table 1 for the full list of surface
physiographic fields; for more details, see IFS documentation [33].

####  2.2.2 ERA5

Climate reanalyses combine observations and modelling to provide calculated
values of a range of climactic variables over time. ERA5 is the fifth-
generation reanalysis from ECMWF. It is produced via 4D-Var data assimilation
of the IFS cycle 41R2 coupled to a land surface model (ECLand, [4]), which
includes lake parametrization by Flake [34] and an ocean wave model (WAM). The
resulting data product provides hourly values of climatic variables across the
atmosphere, land, and ocean at a resolution of approximately 31km with 137
vertical sigma levels up to a height of 80km. Additionally, ERA5 provides
associated uncertainties of the variables at a reduced 63km resolution via a
10-member ensemble of data assimilations. In this work, ERA5 hourly surface
fields at approx. 31km resolution on the cubic octahedral reduced Gaussian
grid (i.e. Tco399) are used. The Gaussian grid’s spacing between latitude
lines is not regular, but lines are symmetrical along the Equator; the number
of points along each latitude line defines longitude lines, which start at
longitude 0 and are equally spaced along the latitude line. In a reduced
Gaussian grid, the number of points on each latitude line is chosen so that
the local east–west grid length remains approximately constant for all
latitudes (here, the Gaussian grid is N320, where N is the number of latitude
lines between a pole and the Equator).

Table 1: Input and target features to all emulators from the data sources. The left column shows the observation-derived static physiographic fields, the middle column ERA5 dynamic physiographic and meteorological fields and the rightmost column ECLand generated dynamic prognostic state variables. Climate fields |  Units |  Atmospheric forcing |  Units |  Prognostic states |  Units
---|---|---|---|---|---
Vegetation cover (low, high) |  |  Total precipitation fraction (convective, stratiform) |  |  Soil water volume (Layers 1-3) |  m3m-3
Type of vegetation (low, high) |  |  Downward radiation (long, short) |  W/m2 |  Soil temperature (Layers 1-3) |  K
Minimum stomatal resistance (low, high) |  |  Seasonal LAI (high, low) |  |  Snow cover fraction |
Roughness length (low, high) |  |  Wind speed (v, u) |  m/s |  |
Urban cover |  |  Surface pressure |  Pa |  |
Lake cover |  |  Skin temperature |  K |  |
Lake depth |  |  |  |  |
Orography (std, filtered) |  m2/s2 |  Specific humidity |  kg/kg |  |
Photosynthesis pathways |  |  Rainfall rate (total) |  kg/m2s |  |
Soil type |  |  Snowfall rate (total) |  kg/m2s |  |
Glacier mask |  |  |  |  |
Permanent wilting point |  |  |  |  |
Field capacity |  |  |  |  |
Cell area |  |  |  |  |

###  2.3 Emulators

We compare the utility of a long-short term memory neural network (LSTM), that
of extreme gradient boosting regression trees (XGB) and that of a feedforward
neural network (that we here refer to as multilayer perceptron, MLP). To
motivate this setup and pave the way for discussing effects of
(hyper-)parameter choices, a short overview of all approaches is given. All
analyses were conducted in Python. XGB was developed in dmlc’s XGBoost python
package111https://xgboost.readthedocs.io/en/stable/python/index.html. The MLP
and LSTM were developed in the PyTorch lightning framework for deep
learning222https://lightning.ai/docs/pytorch/stable/. Neural networks were
trained with the Adam algorithm for stochastic optimization [35]. Model
architectures and algorithmic hyperparameters were selected through Bayesian
hyperparameter optimization with the Optuna framework [36]. The Bayesian
optimization minimizes the neural network validation accuracy, specified here
as mean absolute error (MAE), over a predefined search space for free
hyperparameters with the Tree-structured Parzen Estimator [37]. The resulting
hyperparameter and architecture choices which were used for the different
approaches are listed in the Appendix.

####  2.3.1 MLP

For creation of the MLP emulator we work with a feed-forward neural network
architecture of connected hidden layers with ReLU activations and dropout
layers, model components which are given in detail in the Appendix or in [38].
The MLP was trained with a learning rate scheduler. L2-regularization was
added to the training objective via weight decay. Sizes and width of hidden
layers as well as hyperparameters were selected together in the hyperparameter
optimization procedure. Instead of forecasting absolute prognostic state
variables $z_{t}$, the MLP predicts the 6-hourly increment,
$\frac{\widehat{dz}}{dt}$. It is trained on a stepwise rollout prediction of
future state variables at a pre-defined lead time at given forcing conditions,
see details in the section on optimization.

####  2.3.2 LSTM

LSTMs are recurrent networks that consider long-term dependencies in time
series through gated units with input and forget mechanisms [39]. In
explicitly providing time-varying forcing and state variables, LSTM cell
states serve as long-term memory while LSTM hidden states are the cells’
output and pass on stepwise short-term representations stepwise. In short
notation [20], a one-step ahead forward pass followed by a linear
transformation can be formulated as

$\displaystyle h_{t},\ c_{t}$ $\displaystyle=f\left(x_{t},\ h_{t-1},\
c_{t-1},\ \theta\right)\ $ (1) $\displaystyle{\hat{z}}_{t}$ $\displaystyle=\
{Ah}_{t}+b$

where $h_{t-1}$ denotes the hidden state, i.e. output estimates from the
previous time step, $c_{t-1}$ the cell state from the previous time step, and
$\theta$ the time-invariant model weights. We stacked multiple LSTM cells to
an encoder-decoder model with transfer layers for hidden and cell state
initialization and for transfer to the context vector (see figure 1) [25]. A
lookback $l$ of the previous static and dynamic feature states are passed
sequentially to the first LSTM cells in the encoder layer, while the l
prognostic state variables $z$ initialize the hidden state $h_{0}$ after a
linear embedding. The output of the first LSTM layer cells become the input to
the deeper LSTM layer cells and the last hidden state estimates are the final
output from the encoder. Followed by a non-linear transformation with
hyperbolic tangent activation, the hidden cell states are transformed into a
weighted context vector $s$. Together with the encoder the cell state
($c_{t}$, $s$) initializes the hidden and cell states of the decoder. The
decoder LSTM cells take as input again static and dynamic features
sequentially at lead times $t=1,\ \ldots,\tau$, but not the prognostic states
variables. These are estimated from the sequential hidden states of the last
LSTM layer cells, transformed to target size with a linear forecast head
before prediction. LSTM predicts absolute state variables $z_{t}$ while being
optimized on $z_{t}$ and ${d\hat{z}}_{t}$ simultaneously, see section on
optimization.

Figure 1: LSTM architecture. Blue shaded area indicates the encoder part,
where the model is driven by a lookback $\lambda$ of meteorological forcing
and state variables. The light-blue shaded area indicates the decoder part
that is initialized from the encoding to unroll LSTM forecasts from the
initial time step $t$ up to a flexibly long lead time of $\tau$.

###  2.4 Optimization

####  2.4.1 Loss functions

The basis of the loss function $\mathcal{L}$ for the neural network
optimization was PyTorch’s SmoothL1Loss , a robust loss function that combines
L1-norm and L2-norm and is less sensitive to outliers than pure L1-norm [40].
Based on a pre-defined threshold parameter $\beta$, smooth L1 transitions from
L2-norm to L1-norm above the threshold. SmoothL1Loss $\mathcal{L}$ is defined
as

$\displaystyle\mathcal{L}(\hat{z},\ \ z)$ $\displaystyle=\ 0.5\left(\hat{z}-\
z\right)^{2}\frac{1}{\beta}if|\hat{z}-\ z|<\ \beta\ \text{and}$ (2)
$\displaystyle\mathcal{L}(\hat{z},\ \ z)$ $\displaystyle=\ |\hat{z}-\ z|-\
0.5\ \beta\text{otherwise},$

here with $\beta=1$. All models were trained to minimize the incremental loss
$\mathcal{L}_{s}$ that is the differences between the estimates of the seven
prognostic states increments $\widehat{dz}_{t}$ and the full model’s
prognostic states increments ${dz}_{t}$simultaneously as the sum of losses
over all states. We opted for an loss function equally weighted by variables
to share inductive biases among the non-independent prognostic states [41].
When aggregating over all training lead times $t=1,\ \ldots,\ \tau$,
$\mathcal{L}_{s}$ and grid cells $i=1,\ldots,p$ is

$\displaystyle\mathcal{L}_{s}(\ \widehat{dz},\ dz)=\
\sum_{t}^{\tau}\sum_{i}^{p}{\mathcal{L}_{t}(\ {\widehat{dz}}_{t,i},\
{dz}_{t,i})},$ (3)

Whereas when computing a rollout loss $\mathcal{L}_{r}$ stepwise,

$\displaystyle\mathcal{L}_{r}(\ \widehat{dz},\ z)=\ \frac{1}{\tau\
}\sum_{t}^{\tau}\sum_{i}^{p}{\mathcal{L}_{t}(z_{t-1,i}+\
{\widehat{dz}}_{t,i},\ z_{t,i})}$ (4)

Prognostic state increments are essentially the first differences from one to
the next timestep that are normalized again by the global standard deviation
of the model’s states increments, $s_{\ {}_{dz}}$ before computation of the
loss [keisler_forecasting_2022]. Due to the forecast models’ structural
differences, loss functions were individually adapted.

##### MLP

The combined loss function for the MLP is the sum of the incremental loss
$\mathcal{L}_{s}$ and the rollout loss $\mathcal{L}_{r}$. For the rollout loss
$\mathcal{L}_{r}$, $\mathcal{L}$ was aggregated over grid cells p and
accumulated after an auto-regressive rollout over lead times $\tau$, before
being averaged out by division by $\tau$ [42].

##### LSTM

The combined loss function for the LSTM is the sum of the incremental loss
$\mathcal{L}_{s}$, where the ${d\hat{z}}_{t}$ were derived from
${\hat{z}}_{t}$ after the forward pass, and the loss $\mathcal{L}$ computed on
decoder estimates of prognostic states variables, a functionality that
leverages the potential of our LSTM structure.

##### XGB

Trained only from one to the next time step, i.e. at a lead time of $\tau=1$,
the incremental loss $\mathcal{L}_{s}=\ \mathcal{L}_{r}$. Without a
SmoothL1Loss implementation provided in dmlc’s XGBoost, we trained XGB with
both the Huber-Loss and the default L2-loss. The latter initially providing
better results, we chose the default L2-norm as loss function for XGB with the
regularization parameter $\lambda=1$.

####  2.4.2 Normalization

As prognostic target variables are all lower bounded by zero, we tested both
z-scoring and max-scoring. The latter yielded no significant improvement, thus
we show our results with z-scored target variables. For neural network
training but not for fitting XGB, static, dynamic and prognostic state
variables were all normalized with z-scoring towards the continental or global
mean $\bar{z}$ and unit standard deviation $s_{z}$ as

$z_{t,n}=\ \frac{z_{t,n}-\ \bar{z}}{s_{z}}.$ (5)

Prognostic target state increments were normalized again by the global
standard deviation of increments computing the loss (see section 2.5.1) to
smooth magnitudes of increments [42]. State variables were backtransformed to
original scale before evaluation.

####  2.4.3 Spatial and temporal sampling

Sequences were sampled randomly from the training data set, while validation
happened sequentially. MLP and XGB were trained on all grid cells
simultaneously in both the continental and global setting, while LSTM was
trained on the full continental data set but was limited by GPU memory in the
global task. We overcame this limitation by randomly subsetting grid cells in
the training data into largest possible, equally sized subsets which were then
loaded along with the temporal sequences during the batch sampling.

####  2.4.4 Evaluation

Three scores are used for model validation during the model development phase
and in validating architecture and hyperparameter selection, being the root
mean squared error (RMSE), the mean absolute error (MAE) and the anomaly
correlation coefficient (ACC). First, scores were assessed objectively in
quantifying forecast accuracy of the emulators against ECLand simulations
directly with RMSE and MAE. Doing so, scores were aggregated over lead times,
grid cells or both. The total RMSE was computed as

$\text{RMSE}=\ \sqrt{\frac{\sum_{\tau,\ p}{(z-\ \hat{z})}^{2}}{n}},$ (6)

As the mean absolute error in prognostic state variable prediction over the
total of n grid cells p times lead times $\tau$. Equivalently, MAE was
computed as

$\text{MAE}=\ \frac{\sum_{t,p}{|z-\ \hat{z}|}}{n},$ (7)

Beyond accuracy, the forecast skill of emulators was assessed using
climatology as a benchmark models. The ACC (see below) as index of the long-
term naïve climatology c, i.e. the 6-hourly mean of prognostic state variables
over the last 10 years preceding the test year, i.e. the years 2010 to 2020.
While climatology is a hard-to-beat benchmark specifically in long-term
forecasting, the persistence is a benchmark for short-term forecasting
[pappenberger_how_2015]. For verification against climatology, we compute the
anomaly correlation coefficient (ACC) over lead times as

$\text{ACC}(t)=\frac{\overline{(\hat{z}-c)(z-c)}}{\sqrt{\overline{(\hat{z}-c)^{2}}\,\overline{(z-c)^{2}}}}$
(8)

at each $t=1,…,\tau$ where the overbar denotes averaging over grid cells
$p=i,\ \ldots,\ n$. The nominator now indicates the mean squared skill error
towards climatology and the denominator its variability. ACC is bounded
between 1 and -1, and an ACC of 1 indicates perfect representation of forecast
error variability, an ACC of 0.5 indicates a similar forecast error to that of
the climatology, an ACC of 0 indicates that forecast error variability
dominates and the forecast has no value and an ACC approaching -1 indicates
that the forecast has been very unreliable [43]. ACC is undefined when the
denominator is zero. This is the case either when mean squared emulator or
ECLand anomaly, or both are zero because forecast and climatology perfectly
align, or because they cancel out at summation to the mean.

####  2.4.5 Forecast horizons

Forecast horizons of the emulators are defined by the decomposition of the
RMSE [44] into the emulator’s variability around climatology (i.e. anomaly),
ECLand’s variability around climatology and the covariance of both. The
horizon is the point in time at which the forecast error reaches saturation
level, that is when the covariance of emulator and ECLand anomalies approaches
zero, as does the ACC. We analysed predictive ability and predictability by
computing the ACC for all lead times from 6 hours to approx. one year, i.e.
lead times $t=1,\ldots,\ \tau$, $\tau$ being 1350. As this confounds the
seasonality with the lead time, we compute these for every starting point of
the prediction, requiring two test years (2021 and 2022). Forecast horizons
based on the emulators’ skill in standardized anomaly towards persistence were
equivalently computed but with persistence as a benchmark for shorter time
scales, this was only done for three months, from January to March 2021. The
analysis was conducted on two exemplary regions in northern and southern
Europe that represent very different conditions orography and in prognostic
land surface states, specifically in snow cover. For details on the regions
and on the horizons computed with standardized anomaly skill, see Appendices
A1 and A4 respectively.

##  3 Results

###  3.1 Aggregated performances

##### Europe

All emulators approximated the numerical LSM with high average total
accuracies (all RMSEs < 1.58 and MAEs < 0.84) and confident correlations (all
ACC > 0.72) (see table 2 and figure 2). The LSTM emulator achieved the best
results across all total average scores on the European scale. It decreased
the total average MAE by $\sim$25% towards XGB and by $\sim$37% towards the
MLP and the total average RMSE by $\sim$42% towards XGB and $\sim$38% towards
the MLP. In total average ACC, the LSTM scored 20% higher than the MLP and 15%
than XGB, also being the only emulator that achieved an ACC > 0.9. While the
MLP outperforms XGB in total average RMSE by $\sim$5%, XGB scores better than
the MLP in MAE by $\sim$27%. At variable level, results differentiate into
model specific strengths. In soil water volume, XGB outperforms the neural
network emulators by up to 60% in the first and second layer MAEs towards the
LSTM and up to over 40% for towards the MLP (see table 3). While the
representation of anomalies by specifically the LSTM decreases towards lower
soil layers with an ACC of only 0.6214 at the third soil layer, it remains
consistently higher for XGB with an ACC still > 0.789 at soil layer three.

In soil temperature approximation, LSTM achieves best accuracies at higher
soil levels with up to 7% improvement in MAE towards XGB and ACCs > 0.92, but
XGB outperforms LSTM at the third soil level with a close to 50% improvement
(see table 4). The MLP doesn’t stand out by high scores on the continental
scale. However, in terms of accuracy we found an inverse ranking in the model
development procedure during which LSTM outscored XGB in soil water volume but
struggled with soil temperature approximations, for the interested reader we
refer to the supplementary information. In snow cover approximation, the LSTM
emulator enhances accuracies by over $\sim$50% in MAE towards both the XGB and
the MLP emulator and scores highest in anomaly representation with an ACC of
$\sim$0.87 compared to an ACC of $\sim$0.66 for the MLP and only $\sim$0.74
for the XGB (see table 5).

##### Globe

Score ranking on the global scale varies strongly from the continental scale
(see table 2). In total average accuracies, the MLP outperforms XGB by over
30% and LSTM by up to 25% in RMSE and improves MAE more than 15% towards both.
In anomaly correlation however it scores last, whereas XGB achieves the
highest total average of over 0.75. Consistent with scores on the continental
scale is XGBs high performance in soil temperature (see table 3). It
significantly outperforms the LSTM by $\sim$60% in RMSE and nearly up to 75%
in MAE in all layers and the MLP by up to 50% in MAE at the top layer. Anomaly
persistence for all models degrade visibly towards the lower soil layers,
while that of the LSTM most relative to MLP and XGB. Similar to the
continental scale, XGB also outperforms the other candidates in soil
temperature forecasts in all but the medium layer, where the MLP gets higher
scores in MAE and RMSE but not in ACC (see table 4). LSTM doesn’t stand out
with any scores on the global scale.

Table 2: Emulator total average scores, aggregated over variables, time and space from the European and Global model testings Variable | Model | RMSE | MAE | ACC
---|---|---|---|---
Europe | Globe | Europe | Globe | Europe | Globe
All variables | XGB | 1.575 | 2.611 | 0.695 | 1.601 | 0.765 | 0.755
MLP | 1.486 | 1.699 | 0.832 | 1.189 | 0.728 | 0.569
LSTM | 0.918 | 2.252 | 0.526 | 1.787 | 0.925 | 0.647
Table 3: Emulator average scores on soil water volume forecasts for the European subset, aggregated over space and time from the European and Global model testing. Variable | Layer | Model | RMSE |  | MAE |  | ACC
---|---|---|---|---|---|---|---
|  |  | Europe | Globe |  | Europe | Globe |  | Europe | Globe
Soil water volume | 1 | XGB | 0.013 | 0.015 |  | 0.010 | 0.010 |  | 0.908 | 0.920
|  | MLP | 0.019 | 0.029 |  | 0.015 | 0.023 |  | 0.856 | 0.791
|  | LSTM | 0.029 | 0.048 |  | 0.023 | 0.040 |  | 0.847 | 0.729
| 2 | XGB | 0.011 | 0.012 |  | 0.008 | 0.009 |  | 0.901 | 0.884
|  | MLP | 0.019 | 0.023 |  | 0.014 | 0.018 |  | 0.789 | 0.770
|  | LSTM | 0.029 | 0.050 |  | 0.023 | 0.042 |  | 0.790 | 0.617
| 3 | XGB | 0.015 | 0.014 |  | 0.011 | 0.010 |  | 0.789 | 0.777
| MLP | 0.020 | 0.020 |  | 0.017 | 0.016 |  | 0.576 | 0.667
| LSTM | 0.033 | 0.051 |  | 0.027 | 0.043 |  | 0.621 | 0.475
Table 4: Emulator average scores on soil temperature forecasts for the European subset, aggregated over space and time from the European and Global model testing. Variable | Layer | Model | RMSE |  | MAE |  | ACC
---|---|---|---|---|---|---|---
|  |  | Europe | Globe |  | Europe | Globe |  | Europe | Globe
Soil temperature | 1 | XGB | 1.154 | 4.539 |  | 0.744 | 3.278 |  | 0.806 | 0.769
|  | MLP | 1.628 | 2.606 |  | 1.188 | 2.072 |  | 0.674 | 0.581
|  | LSTM | 0.931 | 3.152 |  | 0.682 | 2.626 |  | 0.938 | 0.735
| 2 | XGB | 0.901 | 2.501 |  | 0.51 | 1.772 |  | 0.812 | 0.797
|  | MLP | 1.134 | 1.851 |  | 0.784 | 1.452 |  | 0.718 | 0.606
|  | LSTM | 0.734 | 2.87 |  | 0.541 | 2.4 |  | 0.928 | 0.699
| 3 | XGB | 0.714 | 1.287 |  | 0.482 | 0.933 |  | 0.722 | 0.711
| MLP | 1.128 | 1.375 |  | 0.821 | 1.071 |  | 0.416 | 0.514
| LSTM | 1.141 | 3.466 |  | 0.918 | 3.002 |  | 0.598 | 0.406
Table 5: Emulator average scores on snow cover forecasts for the European subset, aggregated over space and time from the European and Global model testing. Variable | Layer | Model | RMSE |  | MAE |  | ACC
---|---|---|---|---|---|---|---
|  |  | Europe | Globe |  | Europe | Globe |  | Europe | Globe
Snow Cover | Top | XGB | 8.219 | 9.906 |  | 3.099 | 5.196 |  | 0.746 | 0.707
|  | MLP | 6.449 | 5.995 |  | 2.986 | 3.671 |  | 0.66 | 0.618
|  | LSTM | 3.526 | 6.127 |  | 1.47 | 4.357 |  | 0.877 | 0.698

###  3.2 Spatial and temporal performances

##### Europe

When summarizing temporally aggregated scores as boxplots to a total
distribution over space (see Figure 2, a), the long tails of XGB scores become
visible, whereas the MLP indicates most robustness. This is reflected in the
geographic distribution of scores at the example of ACC (see Figure 2,
bottom), where the area of low anomaly correlation is largest for XGB, ranging
over nearly all northern Scandinavia, while MLP and LSTM have smaller and more
segregated areas of clearly low anomaly correlation. The LSTM shows a
homogenously high ACCs over most of central Europe but the Alps, while also
seems to be challenged in areas of relative to the central Europe extreme
weather conditions at the Norwegian and Spanish coasts.

##### Globe

Similar to the results from the continental analysis, we find again long upper
tails of outliers for XGB in total spatial distribution of accuracies, both in
RMSE and MAE and only few outliers for MLP and LSTM (see Figure 3). The
anomaly correlation distribution changed towards longer lower tails for MLP
and LSTM and a shorter lower tail for XGB. We should, however, take the
results of total average ACC with care as it remains largely undefined in
regions without much noise in snow cover or soil water volume and globally
represents mainly patterns of soil temperature.

Figure 2: a: Total aggregated distributions of (log) scores averaged over lead
times, i.e. displaying the variation among grid cells. b: The distribution of
the anomaly correlation in space on the European subset (b.1: XGB, b.2: MLP,
b.3: LSTM). c: Model forecasts over test year 2021 for grid cell with minimum
and maximum RMSE values (LSTM). Figure 3: a Total average scores, representing
spatial variation among grid cells. b Total average ACC in space. Note that
ACC remained undefined for regions of low signal in snow cover and soil water
volume, see Appendix.

###  3.3 Forecast horizons

Forecast horizons were computed for two European regions, of which the
northern one represents the area of lowest emulators’ skill (see Figure 2,
b.1-3) and the southern one an area of stronger emulators’ skill. Being
strongly correlated with soil water volume, these two regions differ
specifically in their average snow cover fraction (see Figure 4). The
displayed horizons were computed over all prognostic state variables
simultaneously, while their interpretation is related to horizons computed for
prognostic state variables separately, for the figures of which we refer to
the Appendix.

In the North, predictive skill depended on an interaction of how far ahead a
prediction was made (the lead time) and the day of year to which the
prediction was made. In the best case, the LSTM, summer predictions were poor
(light patches in Figure 4 heat maps), but only when initialised in winter.
Or, in other words, one can make good predictions starting in winter, but not
to summer. Vertical structures indicate a systematic model error that appears
at specific initialisation times and that is independent of prediction date,
for example in XGB forecasts that are initialized in May (see Figure 4,
northern region). Diagonal light structures in the heat maps indicate a
temporally consistent error and can be interpreted as physical limits of
system predictability, where the different initial forecast time doesn’t
affect model scores.

All models show stronger limits in predictability and predictive ability in
the northern European region (see Figure 4, left column). MLP and XGB
struggled with representing seasonal variation towards climatology at long
lead times, while LSTM is strongly limited by a systematic error in certain
regions. Initializing the forecast the 1 January 2021, MLP drops below an ACC
of 80% repeatedly from initialization on and then to an ACC below 10% at the
beginning of May. LSTMs performance is more robust in the beginning of the
year but depletes strongly later to less than 10% ACC in mid May. On the one
hand, this represents two different characteristics of model errors: MLP
forecasts for snow cover fraction are less than zero for some grid cells while
LSTM forecasts for snow cover fraction remain falsely at very high levels for
some grid cells, not predicting the snowmelt in May (see Appendix, A4.1). On
the other hand, this represents a characteristic error due to change in
seasonality: the snowmelt in this region in May happens abruptly and all
emulators repeatedly over- or underpredict the exact date.

Figure 4: Emulator forecast skill horizons in two European subregions,
aggregated over prognostic state variables. Scores are computed with the
anomaly correlation coefficient (ACC) at 6-hourly lead times (y-axis) over
approx. one year, displayed as a function of the initial forecast time
(x-axis). As horizon we define the time at which the forecast has no value at
all, i.e. when ACC is 0 (or below 10%). The diagonal dashed lines indicate the
day of the test year 2021 as labelled on the x-axis, the arrows indicate where
forecasts reach the second test year 2022.

##  4 Discussion

In the comparative analysis of emulation approaches for land surface
forecasting, three primary models—LSTM (Long Short-Term Memory networks), MLP
(Multi-Layer Perceptrons), and XGB (Extreme Gradient Boosting)—have been
evaluated to understand their effectiveness across different operational
scenarios. While all emulators achieved high predictive scores, models differ
in their demand of computational resources [45] and each one offers unique
advantages and faces distinct challenges, impacting their suitability for
various forecasting tasks. With this work we want to present the first steps
towards enabling quick offline experimentation on the land surface with
ECMWF’s land surface scheme ECLand and decreasing computational demands, i.e.
in the coupled data assimilation.

###  4.1 Approximation of prognostic land surface states

The total evaluation scores of our emulators indicate good agreement with
ECLand simulations. Among the seven individual prognostic land surface states,
emulators achieve notably different scores and in the transfer from the high-
resolution continental to the low-resolution global scale, their performance
ranking change. On average, neural network performances degrade towards the
deeper soil layers, while XGB scores remain relatively stable. Also, the
neural networks scores drop in the extrapolation from continental to global
scale, while XGB scores also for this task remain constantly high. In a way,
these findings are not surprising. It is known that neural networks are highly
sensitive to selection bias [17] and tuning of hyper-parameters [46],
suboptimal choices of which may destabilise variance in predictive skill.
Previous and systematic comparisons of XGB and deep neural networks have
demonstrated that neural networks can hardly be transferred to new data sets
without performance loss [16]. On tabular data, XGB still outperforms neural
networks in most cases [17], unless these models are strongly regularized
[47]. The disadvantage of neural networks might lay in the rotational
invariance of MLP-like architectures, due to which information about the data
orientation gets lost, as well as in their instability regarding uninformative
input features [17].

Inversely to expectations and preceding experiments, on the European data set
relative to the two other models the LSTM scored better in the upper layer
soil temperatures than in forecasting soil water volume and decreased in
scores towards lower layers with slower processes. For training on
observations, the decreasing LSTM predictive accuracy for soil moisture with
lead time is discussed [48], but reasons arising from the engineering side
remain unclear. In an exemplary case of a single-objective, deterministic
streamflow forecast, a decrease in recurrent neural network performance has
been related with an increasing coefficient of variation [49]. In our European
subregions, the signal-to-noise ratio of the prognostic state variables
(computed as the averaged ratio of mean and standard deviation) is up to ten
times higher in soil temperature than in soil water volume states (see
Appendix, A2.1). While a small signal of the latter may induce instability in
scores, it does not explain the decreasing performance towards deeper soil
layers with slow processes, where we expected an advantage of the long-term
memory.

Stein’s paradox tells us that joint optimization may lead to better results if
the target is multi-objective, but not if we are interested in single targets
[50, 41]. While from a process perspective multi-objective scores are less
meaningful than single ones, this is what we opted for due to efficiency. The
unweighted linear loss combination might be suboptimal in finding effective
parameters across all prognostic state variables [51, 41], yet being strongly
correlated, we deemed their manual weighting inappropriate. An alternative to
this provides adaptive loss weighting with gradient normalisation [51].

###  4.2 Evaluation in time and space

We used aggerated MAE and RMSE accuracies as a first assessment tool to
conduct model comparison, but score aggregation hides model specific spatio-
temporal residual patterns. Further, both scores are variance dependent,
favouring low variability in model forecasts even though this may not be
representative of the system dynamic [52]. Assessing the forecast skill over
time as the relative proximity to a subjectively chosen benchmark helps
disentangling areas of strengths and weaknesses in forecasting with the
emulators [26]. The naïve 6-hourly climatology as benchmark highlights periods
where emulators long-range forecasts on the test year are externally limited
by seasonality, i.e. system predictability, and where they are internally
limited by model error, i.e. the model’s predictive ability. Applying this
strategy in two exemplary European subregions showed that all emulators
struggle most in forecasting the period from late summer to autumn, unless
they are initialized in summer (see Figure 4). Because forecast quality is
most strongly limited by snow cover (see Appendix, A4.1), we interpret this as
the unpredictable start of snow fall in autumn. External predictability
limitations seem to affect the LSTM overall less than the two other models,
and specifically XGB drifts at long lead times. From a geographical
perspective inferred from the continental scale, emulators struggle in
forecasting prognostic state variables in regions with complicated orography
and strong environmental gradients. XGB scores vary seemingly random in space,
while neural networks scores exhibit spatial autocorrelation. A meaningful
inference about this, however, can only be conducted in assessing model
sensitivities to physiographic and meteorological fields through gradients and
partial dependencies. While the goal of this work is to introduce our approach
to emulator development, we envision this for follow-up analyses.

###  4.3 Emulation with memory mechanisms

Without much tuning, XGB challenges both LSTM and MLP for nearly all variables
(see tables 2-4). In training on observations for daily short-term and real-
time rainfall-runoff prediction, XGB and LightXGB were shown before to equally
performed as, or outperformed LSTMs [51, 45]. Nevertheless, models with memory
mechanism such as the encoder-decoder LSTM remain a promising approach for
land surface forecasting regarding their differentiability [18], their
flexible extension of lead times, for exploring the effect of long-term
dependencies or for inference from the context vector that may help
identifying the process relevant climate fields [20].

In our LSTM architecture, we assume that our model is well defined in that the
context vector perfectly informs the hidden decoder states. If that assumption
is violated, potential strategies are to create a skip-connection between
context vector and forecast head, or to consider input of time-lagged
variables or self-attention mechanisms [53]. With attention, the context
vector becomes a weighted sum of alignments that relates neighbouring
positions of a sequence, a feature that could be leveraged for forecasting
quick processes such as snow cover or top-level soil water volume.

Comparing average predictive accuracies across different training lead times
indicates that training at longer lead times may enhance short-term accuracy
of the LSTM at the cost of training runtime (see Appendix, A2). A superficial
exploration of encoder length indicates no visible improvement on target
accuracies if not a positive tendency towards shorter sequences. This needs an
extended analysis for understanding, yet without a significant improvement by
increased sequence length, GRU cells might provide a simplified and less
parameterized alternative to LSTM cells. They were found to perform equally
well on streamflow forecast performance before, while reaching higher
operational speed [49].

###  4.4 Emulators in application

LSTM networks with a decoder structure are valued for their flexible and fast
lead time evaluation, which is crucial in applications where forecast
intervals are not consistent. The structure of LSTM is well-suited for
handling sequential data, allowing it to perform effectively over different
temporal scales [39]. They provide access to gradients, which facilitates
inference, optimization and usage for coupled data assimilation [18, 54].
Nevertheless, the complexity of LSTMs introduces disadvantages: Despite their
high evaluation speed and accuracy under certain conditions, they require
significant computational resources and long training times. They are also
highly sensitive to hyperparameters, making them challenging to tune and slow
to train, especially with large datasets.

MLP models stand out for their implementation, training and evaluation speed
with yet rewarding accuracy, making them a favourable choice for scenarios
that require rapid model deployment. They are tractable and easy to handle,
with a straightforward setup that is less demanding computationally than more
complex models. MLPs also allow for access to gradients, aiding in incremental
improvements during training and quick inference [18]. Despite these
advantages, MLPs face challenges with memory scaling during training at fixed
lead times, which can hinder their applicability in large-scale or high-
resolution forecasting tasks.

XGB models are highly regarded for their robust performance with minimal
tuning, achieving high accuracy not only in sample applications, but also in
transfer to unseen problems [16, 17]. Their simplicity makes them easy to
handle, even for users with limited technical expertise in machine learning.
However, the slow evaluation speed of XGB becomes apparent as dataset
complexity and size increase. Although generally more interpretable than deep
machine learning tools, XGB is not differentiable, limiting its application in
coupled data assimilation [18, 54] even though research on differentiable
trees is ongoing [55].

##  5 Conclusion

In conclusion, the choice between LSTM, MLP, and XGB models for land surface
forecasting depends largely on the specific requirements of the application,
including the need for speed, accuracy, and ease of use. Each model’s
computational demands, flexibility, and operational overhead must be carefully
considered to optimize performance and applicability in diverse forecasting
environments. When it comes to accuracy, combined model ensembles of XGB and
neural networks have been shown to yield the best results [16], but accuracy
alone will not determine a single best approach [46]. Our comparative
assessment underscores the importance of selecting the appropriate emulation
approach based on a clear understanding of each model’s strengths and
limitations in relation to the forecasting tasks at hand. By developing the
emulators for ECMWF’s numerical land surface scheme ECLand, we path the way
towards a physics-informed ML-based land surface model that on the long run
can be parametrized with observations. We also provide a pretrained model
suite to improve land surface forecasts and future land reanalyses.

## Code and data availability

Code for this analysis can be found here: https://github.com/MWesselkamp/land-
surface-emulation. Data is available on request.

## Author contribution

MW, MCha, EP, FP and GB conceived the study. MW and EP conducted the analysis.
MW, MCha, MK, EP discussed and took technical decisions. SB advised on process
decisions. MW, MCho and FP wrote the manuscript. EP, CFD, FP reviewed the
analysis and/or manuscript.

## Competing interest

The authors declare that they have no conflict of interest.

## Acknowledgements

This work profited from discussion with Linus Magnusson, Patricia de Rosnay,
Sina R. K. Farhadi and Karan Ruparell and many more. MW thankfully
acknowledges ECMWF for providing two research visit stipendiates over the
course of the collaboration. EP was funded by the CERISE project (grant
agreement No101082139) funded by the European Union. Views and opinions
expressed are however those of the authors only and do not necessarily reflect
those of the European Union or the Commission. Neither the European Union nor
the granting authority can be held responsible for them. ChatGPT version 4.0
was used for coding support.

## References

  * [1] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia.  Learning skillful medium-range global weather forecasting.  Science, 382(6677):1416–1421, December 2023.
  * [2] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian.  Pangu-Weather: A 3D High-Resolution Model for Fast and Accurate Global Weather Forecast.  2022\.
  * [3] Simon Lang, Mihai Alexe, Matthew Chantry, Jesper Dramsch, Florian Pinault, Baudouin Raoult, Mariana C. A. Clare, Christian Lessig, Michael Maier-Gerber, Linus Magnusson, Zied Ben Bouallègue, Ana Prieto Nemesio, Peter D. Dueben, Andrew Brown, Florian Pappenberger, and Florence Rabier.  AIFS - ECMWF’s data-driven forecasting system, June 2024.
  * [4] Souhail Boussetta, Gianpaolo Balsamo, Gabriele Arduini, Emanuel Dutra, Joe McNorton, Margarita Choulga, Anna Agustí-Panareda, Anton Beljaars, Nils Wedi, Joaquín Munõz-Sabater, Patricia De Rosnay, Irina Sandu, Ioan Hadade, Glenn Carver, Cinzia Mazzetti, Christel Prudhomme, Dai Yamazaki, and Ervin Zsoter.  ECLand: The ECMWF land surface modelling system.  Atmosphere, 12(6):723, June 2021.
  * [5] Patricia De Rosnay, Gianpaolo Balsamo, Clément Albergel, Joaquín Muñoz-Sabater, and Lars Isaksen.  Initialisation of land surface variables for numerical weather prediction.  Surveys in Geophysics, 35(3):607–621, May 2014.
  * [6] Joaquín Muñoz-Sabater, Emanuel Dutra, Anna Agustí-Panareda, Clément Albergel, Gabriele Arduini, Gianpaolo Balsamo, Souhail Boussetta, Margarita Choulga, Shaun Harrigan, Hans Hersbach, Brecht Martens, Diego G. Miralles, María Piles, Nemesio J. Rodríguez-Fernández, Ervin Zsoter, Carlo Buontempo, and Jean-Noël Thépaut.  ERA5-Land: A state-of-the-art global reanalysis dataset for land applications.  Earth System Science Data, 13(9):4349–4383, September 2021.
  * [7] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, and Prabhat.  Deep learning and process understanding for data-driven Earth system science.  Nature, 566(7743):195–204, February 2019.
  * [8] David Machac, Peter Reichert, and Carlo Albert.  Emulation of dynamic simulators with application to hydrology.  Journal of Computational Physics, 313:352–366, May 2016.
  * [9] Istem Fer, Ryan Kelly, Paul R. Moorcroft, Andrew D. Richardson, Elizabeth M. Cowdery, and Michael C. Dietze.  Linking big models to big data: Efficient ecosystem model calibration through Bayesian model emulation.  Biogeosciences, 15(19):5801–5830, October 2018.
  * [10] Matthew Chantry, Sam Hatfield, Peter Duben, Inna Polichtchouk, and Tim Palmer.  Machine learning emulation of gravity wave drag in numerical weather forecasting.  2021\.
  * [11] Lizao Li, Rob Carver, Ignacio Lopez-Gomez, Fei Sha, and John Anderson.  SEEDS: Emulation of Weather Forecast Ensembles with Diffusion Models.  2023\.
  * [12] Marieke Wesselkamp, Niklas Moser, Maria Kalweit, Joschka Boedecker, and Carsten F. Dormann.  Process-guidance improves predictive performance of neural networks for carbon turnover in ecosystems.  2022\.
  * [13] Shruti Nath, Quentin Lejeune, Lea Beusch, Sonia I. Seneviratne, and Carl-Friedrich Schleussner.  MESMER-M: An Earth system model emulator for spatially resolved monthly temperature.  Earth System Dynamics, 13(2):851–877, April 2022.
  * [14] Evan Baker, Anna B. Harper, Daniel Williamson, and Peter Challenor.  Emulation of high-resolution land surface models using sparse Gaussian processes with application to JULES.  Geoscientific Model Development, 15(5):1913–1929, March 2022.
  * [15] David Meyer, Sue Grimmond, Peter Dueben, Robin Hogan, and Maarten Van Reeuwijk.  Machine Learning Emulation of Urban Land Surface Processes.  Journal of Advances in Modeling Earth Systems, 14(3):e2021MS002744, March 2022.
  * [16] Ravid Shwartz-Ziv and Amitai Armon.  Tabular Data: Deep Learning is Not All You Need.  2021\.
  * [17] Léo Grinsztajn, Edouard Oyallon, and Gaël Varoquaux.  Why do tree-based models still outperform deep learning on tabular data?, 2022.
  * [18] Sam Hatfield, Matthew Chantry, Peter Dueben, Philippe Lopez, Alan Geer, and Tim Palmer.  Building Tangent-Linear and Adjoint Models for Data Assimilation With Neural Networks.  Journal of Advances in Modeling Earth Systems, 13(9):e2021MS002521, September 2021.
  * [19] Frederik Kratzert, Daniel Klotz, Guy Shalev, Günter Klambauer, Sepp Hochreiter, and Grey Nearing.  Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets.  Hydrology and Earth System Sciences, 23(12):5089–5110, December 2019.
  * [20] Thomas Lees, Steven Reece, Frederik Kratzert, Daniel Klotz, Martin Gauch, Jens De Bruijn, Reetik Kumar Sahu, Peter Greve, Louise Slater, and Simon J. Dadson.  Hydrological concept formation inside long short-term memory (LSTM) networks.  Hydrology and Earth System Sciences, 26(12):3079–3101, June 2022.
  * [21] Jacob A. Zwart, Samantha K. Oliver, William David Watkins, Jeffrey M. Sadler, Alison P. Appling, Hayley R. Corson-Dosch, Xiaowei Jia, Vipin Kumar, and Jordan S. Read.  Near-term forecasts of stream temperature using deep learning and data assimilation in support of management decisions.  JAWRA Journal of the American Water Resources Association, 59(2):317–337, April 2023.
  * [22] Alberto Bassi, Marvin Höge, Antonietta Mira, Fabrizio Fenicia, and Carlo Albert.  Learning Landscape Features from Streamflow with Autoencoders, February 2024.
  * [23] Frederik Kratzert, Mathew Herrnegger, Daniel Klotz, Sepp Hochreiter, and Günter Klambauer.  NeuralHydrology – Interpreting LSTMs in Hydrology.  2019\.
  * [24] Peter Van Katwyk, Baylor Fox-Kemper, Hélène Seroussi, Sophie Nowicki, and Karianne J. Bergen.  A Variational LSTM Emulator of Sea Level Contribution From the Antarctic Ice Sheet.  Journal of Advances in Modeling Earth Systems, 15(12):e2023MS003899, December 2023.
  * [25] Grey Nearing, Deborah Cohen, Vusumuzi Dube, Martin Gauch, Oren Gilon, Shaun Harrigan, Avinatan Hassidim, Daniel Klotz, Frederik Kratzert, Asher Metzger, Sella Nevo, Florian Pappenberger, Christel Prudhomme, Guy Shalev, Shlomo Shenzis, Tadele Yednkachw Tekalign, Dana Weitzner, and Yossi Matias.  Global prediction of extreme floods in ungauged watersheds.  Nature, 627(8004):559–563, March 2024.
  * [26] F. Pappenberger, M.H. Ramos, H.L. Cloke, F. Wetterhall, L. Alfieri, K. Bogner, A. Mueller, and P. Salamon.  How do I know if my forecasts are better? Using benchmarks in hydrological ensemble prediction.  Journal of Hydrology, 522:697–713, March 2015.
  * [27] Gianpaolo Balsamo, Souhail Boussetta, Emanuel Dutra, Anton Beljaars, and Pedro Viterbo.  Evolution of land-surface processes in the IFS.  (127), 2011.
  * [28] ECMWF.  IFS Documentation CY43R3 - Part IV: Physical processes.  2017\.
  * [29] Tom Kimpson, Margarita Choulga, Matthew Chantry, Gianpaolo Balsamo, Souhail Boussetta, Peter Dueben, and Tim Palmer.  Deep learning for quality control of surface physiographic fields using satellite Earth observations.  Hydrology and Earth System Sciences, 27(24):4661–4685, December 2023.
  * [30] Pedro Viterbo.  Land_surface_processes.  2002\.
  * [31] Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, Adrian Simmons, Cornel Soci, Saleh Abdalla, Xavier Abellan, Gianpaolo Balsamo, Peter Bechtold, Gionata Biavati, Jean Bidlot, Massimo Bonavita, Giovanna De Chiara, Per Dahlgren, Dick Dee, Michail Diamantakis, Rossana Dragani, Johannes Flemming, Richard Forbes, Manuel Fuentes, Alan Geer, Leo Haimberger, Sean Healy, Robin J. Hogan, Elías Hólm, Marta Janisková, Sarah Keeley, Patrick Laloyaux, Philippe Lopez, Cristina Lupu, Gabor Radnoti, Patricia De Rosnay, Iryna Rozum, Freja Vamborg, Sebastien Villaume, and Jean-Noël Thépaut.  The ERA5 global reanalysis.  Quarterly Journal of the Royal Meteorological Society, 146(730):1999–2049, July 2020.
  * [32] Margarita Choulga, Ekaterina Kourzeneva, Gianpaolo Balsamo, Souhail Boussetta, and Nils Wedi.  Upgraded global mapping information for earth system modelling: An application to surface water depth at the ECMWF.  Hydrology and Earth System Sciences, 23(10):4051–4076, October 2019.
  * [33] ECMWF.  IFS Documentation CY48R1 - Part IV: Physical Processes.  2023\.
  * [34] Dmitrii Mironov and Jürgen Helmert.  Parameterization of Lakes in NWP and Climate Models.
  * [35] Diederik P. Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization, 2017.
  * [36] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A Next-generation Hyperparameter Optimization Framework.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2623–2631, Anchorage AK USA, July 2019. ACM.
  * [37] Yoshihiko Ozaki, Yuki Tanigaki, Shuhei Watanabe, Masahiro Nomura, and Masaki Onishi.  Multiobjective Tree-Structured Parzen Estimator.  Journal of Artificial Intelligence Research, 73:1209–1250, April 2022.
  * [38] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  Adaptive Computation and Machine Learning. The MIT Press, Cambridge, Massachusetts, 2016.
  * [39] Sepp Hochreiter and Jürgen Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, November 1997.
  * [40] Ross Girshick.  Fast R-CNN.  2015\.
  * [41] Ozan Sener and Vladlen Koltun.  Multi-Task Learning as Multi-Objective Optimization, 2018.
  * [42] Ryan Keisler.  Forecasting Global Weather with Graph Neural Networks, February 2022.
  * [43] ECMWF.  Forecast User Guide.  Anomaly Correlation Coefficient.
  * [44] Lisa K. Bengtsson, Linus Magnusson, and Erland Källén.  Independent Estimations of the Asymptotic Variability in an Ensemble Forecast System.  Monthly Weather Review, 136(11):4105–4112, November 2008.
  * [45] Zhongjie Cui, Xiaoxia Qing, Hongxiang Chai, Senxiong Yang, Ying Zhu, and Feifei Wang.  Real-time rainfall-runoff prediction using light gradient boosting machine coupled with singular spectrum analysis.  Journal of Hydrology, 603:127124, December 2021.
  * [46] Xavier Bouthillier, Pierre Delaunay, Mirko Bronzi, Assya Trofimov, Brennan Nichyporuk, Justin Szeto, Naz Sepah, Edward Raff, Kanika Madan, Vikram Voleti, Samira Ebrahimi Kahou, Vincent Michalski, Dmitriy Serdyuk, Tal Arbel, Chris Pal, Gaël Varoquaux, and Pascal Vincent.  Accounting for Variance in Machine Learning Benchmarks, March 2021.
  * [47] Arlind Kadra, Marius Lindauer, Frank Hutter, and Josif Grabocka.  Well-tuned Simple Nets Excel on Tabular Datasets, November 2021.
  * [48] Pingki Datta and Salah A. Faroughi.  A multihead LSTM technique for prognostic prediction of soil moisture.  Geoderma, 433:116452, May 2023.
  * [49] Yuxue Guo, Xinting Yu, Yue-Ping Xu, Hao Chen, Haiting Gu, and Jingkai Xie.  AI-based techniques for multi-step streamflow forecasts: Application for multi-objective reservoir operation optimization and performance assessment.  Hydrol. Earth Syst. Sci., 2021.
  * [50] W. James and Charles Stein.  Estimation with Quadratic Loss.  In Samuel Kotz and Norman L. Johnson, editors, Breakthroughs in Statistics, pages 443–460. Springer New York, New York, NY, 1992.
  * [51] Zhao Chen, Vijay Badrinarayanan, Chen-Yu Lee, and Andrew Rabinovich.  GradNorm: Gradient Normalization for Adaptive Loss Balancing in Deep Multitask Networks.  2017\.
  * [52] A. Thorpe, Peter Bauer, L. Magnusson, and David Richardson.  An evaluation of recent performance of ECMWF’s forecasts.  2013\.
  * [53] Xi Chen, Jiaxu Huang, Zhen Han, Hongkai Gao, Min Liu, Zhiqiang Li, Xiaoping Liu, Qingli Li, Honggang Qi, and Yonggui Huang.  The importance of short lag-time in the runoff forecasting model based on long short-term memory.  Journal of Hydrology, 589:125359, October 2020.
  * [54] Patricia De Rosnay, Philip Browne, Eric De Boisséson, David Fairbairn, Yoichi Hirahara, Kenta Ochi, Dinand Schepers, Peter Weston, Hao Zuo, Magdalena Alonso-Balmaseda, Gianpaolo Balsamo, Massimo Bonavita, Niels Borman, Andy Brown, Marcin Chrust, Mohamed Dahoui, Giovanna Chiara, Stephen English, Alan Geer, Sean Healy, Hans Hersbach, Patrick Laloyaux, Linus Magnusson, Sébastien Massart, Anthony McNally, Florian Pappenberger, and Florence Rabier.  Coupled data assimilation at ECMWF: Current status, challenges and future developments.  Quarterly Journal of the Royal Meteorological Society, 148(747):2672–2702, July 2022.
  * [55] Sergei Popov, Stanislav Morozov, and Artem Babenko.  Neural Oblivious Decision Ensembles for Deep Learning on Tabular Data, 2019.