Regional data-driven weather modeling with a global stretched-grid
Thomas Nils Nipen Håvard Homleid Haugen Magnus Sikora Ingstad Even Marius Nordhagen Aram Farhad Shafiq Salihi Paulina Tedesco Ivar Ambjørn Seierstad Jørn Kristiansen Simon Lang Mihai Alexe Jesper Dramsch Baudouin Raoult Gert Mertes Matthew Chantry



[
\begin@twocolumnfalse
[
    [
    September 9, 2024
=====================

§ ABSTRACT 

A ddm suitable for regional weather forecasting applications is presented. The model extends the aifs by introducing a stretched-grid architecture that dedicates higher resolution over a regional area of interest and maintains a lower resolution elsewhere on the globe. The model is based on gnns, which naturally affords arbitrary multi-resolution grid configurations.

The model is applied to short-range weather prediction for the Nordics, producing forecasts at 2.5km spatial and 6h temporal resolution. The model is pre-trained on 43 years of global ERA5 data at 31km resolution and is further refined using 3.3 years of 2.5km resolution operational analyses from the meps. The performance of the model is evaluated using surface observations from measurement stations across Norway and is compared to short-range weather forecasts from meps. The ddm outperforms both the control run and the ensemble mean of meps for 2m temperature. The model also produces competitive precipitation and wind speed forecasts, but is shown to underestimate extreme events.






\end@twocolumnfalse
]

[1]Norwegian Meteorological Institute, Oslo, Norway
[2]ECMWF, Reading, UK
[3]ECMWF, Bonn, Germany

^*Equal contribution






§ BACKGROUND

ddms are rapidly emerging as an alternative to nwp models. Global ddms have been shown to produce weather forecasts that are competitive in accuracy with state-of-the-art global nwp models <cit.>. Unlike nwp models that simulate the weather by solving physical equations, ddms are typically based on artificial neural networks with parameters that are optimized by training on historical reanalyses. Although training costs are significant, ddms can be several orders of magnitude cheaper to run than nwp models of similar resolution, due to the efficient mapping of neural networks to gpus and their ability to use a much longer timestep when integrating forward in time.

This progress has been made possible by the availability of free and reliable datasets like ERA5 <cit.>, which provides a long archive of consistent reanalyses at approximately 31km horizontal grid resolution. Several global ddms have been trained using ERA5, including FourCastNet <cit.>, Pangu-Weather <cit.>, GraphCast <cit.>, FuXi <cit.>, and the aifs <cit.>.



The success of global ddms is driving efforts to develop ddms for regional purposes. In their mission to save lives and protect property, nmhses, such as the metno, currently rely on high-resolution regional nwp models to deliver real-time weather forecasts to the public and weather-affected industries.

Additionally, many nmhses produce regional reanalyses, with consistent high-resolution data spanning long time periods. These high-resolution datasets, often at a kilometer scale, offer an enormous potential for training regional ddms. Due to their high resolution, these datasets provide additional information at finer spatial scales than ERA5 can represent.

Regional data-driven modeling efforts have so far focused on driving a lam with initial conditions from a regional analysis and lateral boundaries from an external global model <cit.>. This is analogous to the common separation of regional and global models in nwp. 



This article presents a different approach, where a global model is developed with higher resolution over the region of interest and lower resolution elsewhere. The model is initialized from both global and high-resolution regional analyses. The goal of this stretched-grid approach is for weather systems to move seamlessly from the global domain into the regional domain and out again without any explicit treatment of the boundary between the domains. The model is based on gnns <cit.>, which is also the architecture used in GraphCast and aifs, and was introduced to weather prediction by <cit.>. This flexible architecture allows arbitrary grids to be used, supporting our need to have variable resolution across the globe.

We target the data-driven model for use by metno's public weather forecasting app Yr (<https://www.yr.no>). Given Norway's intricate topography and coastline, the model must be capable of providing high-resolution forecasts. Currently, short-range weather forecasts on Yr for the Nordics are provided by the MetCoOp Ensemble Prediction System (MEPS) <cit.> with a horizontal resolution of 2.5km. The performance of the model is therefore compared to the operational MEPS and the evaluation of the model is done with a focus on the needs of our end-users.    


This article is organized as follows: Section <ref> describes the stretched-grid approach to regional data-driven modeling. Section <ref> describes the hyper-parameters we chose and how the model was trained. Section <ref> evaluates the data-driven model against current operational nwp systems. Section <ref> draws conclusions and presents future work.



§ STRETCHED-GRID MODEL




gnn-based models, such as GraphCast and aifs, use an encoder-processor-decoder architecture (<ref>). Weather parameters are provided to the model on an input grid (<ref>a,b). Typically, this includes data on multiple pressure levels for the current time and one or more recent times in the past. This data is then encoded into a compact latent representation of the state of the atmosphere (<ref>c). The latent space is typically on a lower resolution grid, which we will call the processor mesh (<ref>d). This latent state is advanced forward one time step by the processor. The updated state is then decoded back to the weather parameters on the input grid (<ref>e). This process can be repeated to create forecasts of arbitrary length, which is called rollout.

Graphs are central building blocks of gnns, consisting of a set of nodes connected by edges. Nodes are used to represent the points in the input grid and in the processor mesh. As the nodes in the graphs are represented by a one-dimensional vector (i.e. is not required to be on a regular grid), arbitrary grid configurations can be defined. Graphs are well suited for a problem with variable resolution because the structural information in the model is encoded in the model features while the model weights do not depend explicitly on the shape of the graph. This allows us to dedicate a finer grid spacing in parts of the domain, which we call a stretched grid. Although the stretched-grid approach can be used with any number of different resolution domains, we focus here on applications with a single global and a single regional domain.

Graphs are used by the encoder, processor, and decoder. These components update the nodes they are connected to by processing information from any node connected by the edges. This update is called message passing. To aggregate the information from the contributing nodes, we use a graph transformer, which is a generalization of the transformer architecture <cit.> to graphs, where message passing is done via a multi-head attention mechanism <cit.>, the same GNN architecture used in the encoder and decoder of AIFS <cit.>. Our motivation for choosing the graph transformer over a pure transformer is that the graph transformer explicitly embeds directional and positional information via the edge features into the calculation of attention scores, which can be important for a multi-resolution model. 

The stretched-grid model is built using the Anemoi framework first created by ecmwf and now co-developed by several meteorological organisations across Europe, including metno (<https://anemoi.ecmwf.int/>).



 §.§ Grid and processor mesh design


The input grid is created to match the exact grid points from the two input models. Grid points from the global model that overlap with the regional model are removed, which is a configuration referred to as cutout. To construct the processor mesh, we follow <cit.> and start with an initial regular icosahedron with 20 faces, creating a graph with 12 nodes and 30 edges. This icosahedron is gradually refined by adding nodes and edges that divide each triangular face into four smaller triangles (<ref>f). When a layer of refinement is added to the mesh, we keep the edges of the original lower-resolution mesh so that different levels of refinement facilitate communication on different length scales in the final multi-layer mesh. We add extra levels of refinement over the regional domain, with the aim of keeping the number of edges from the input grid to each mesh node constant between the regional and global domains.

The encoder maps data on the input grid to a hidden representation that encodes the structural information in the graph. The first step is the source node embedding, which consists of a linear layer that transforms the raw features of the source nodes (in the grid) into a latent space representation. Similarly, the destination node embedding transforms the raw features from the destination nodes (in the mesh) to a hidden dimension. Layer normalization is applied to both source and destination node embeddings. The embedded source and destination nodes, along with processed edge attributes, are passed through a graph transformer block to perform the core processing. Each processor mesh node is connected to the nearest 12 input grid points (<ref>c). 



 §.§ Processor

The processor evolves the latent space representation of the data on the processor mesh through a series of message-passing steps. Each message-passing step is done through a graph transformer block where every node receives information from all one-hop neighbors in the multi-layer mesh. Consequently, communication on different mesh resolutions happens simultaneously <cit.> with spatial relations embedded in the multi-head attention mechanism. In total, the processor consists of 16 such message-passing steps, in contrast to one in the encoder and decoder.



 §.§ Decoder

The decoder maps the latent space representation back to weather parameters on the input grid. To achieve this, the embedding process performed by the encoder is reversed by a graph transformer block, followed by embedding of the input features via a skip connection. Each point in the output grid receives data from its three nearest processor mesh nodes (<ref>e).



§ MODEL TRAINING







 §.§ Data

To train the stretched-grid model, we use both global and regional NWP analyses. As global data, we use ERA5 reanalysis <cit.> on its native 31km resolution reduced Gaussian grid for the period 1979–2022. We also use the operational analyses of the ifs <cit.>, interpolated to the ERA5 resolution for the time period 2020–2024, as this is what the model will use as input when run operationally.

For the regional domain, we use the analysis of the meps <cit.> with 2.5km horizontal grid spacing, which covers Norway, Sweden, Finland, Denmark, and the Baltic countries. As the model underwent substantial changes, including a change to the grid configuration, on February 5, 2020, we only use analyses after this time.

<ref> summarizes the input and output variables we used from these datasets.

When performing hyper-parameter exploration, we used any data available in the time period 1979-01-01 to 2022-05-31 for training, reserving 2022-06-01 to 2023-05-31 for validation to help select the best model. We reserved a separate period for independent testing of the best model for 2023-06-01 to 2024-05-31, with results presented in <ref>.



 §.§ Training configuration



Many of the settings used for the model follow <cit.>, including two input time steps, the choice of loss function (squared error), per-variable weights, and optimizer. We used 1024 channels in the encoder, processor, and decoder, resulting in a total of 246 million trainable parameters.

 

Early results showed that training a stretched-grid model using only 2-3 years of data leads to forecasts with poor synoptic developments. Therefore, we pre-trained the model on the 43 years of ERA5 data, before fine-tuning on the combined datasets. A consequence of this is that we cannot use any learnable features (<cit.>), as these are graph-specific.

We used a four-stage training procedure, as illustrated in <ref> and with full details in <ref>. After each stage, the scheduler was reset with a smaller starting learning rate to prevent catastrophic forgetting.

In stage A, the model was pre-trained on ERA5 upsampled to 100 km resolution. This is a relatively cheap procedure that subsequently speeds up the training of Stage B, where we trained on ERA5's full resolution.

In Stage C, we switched to using IFS upsampled to the native ERA5 resolution and added MEPS at 2.5km resolution in a stretched-grid configuration. In <cit.>, each output grid point's contribution to the loss is made proportional to the area that the point covers. We found that increasing the contribution of the grid points in the regional domain improved verification scores for the regional domain. In particular, we let the regional domain contribute 33% to the total loss value despite it covering only 1.2% of the earth's surface.

Finally, in Stage D we perform a 24 h auto-regressive rollout training to improve scores over longer lead times. This is done by incrementally increasing the rollout from two to four 6h time steps with 100 iterations each. Here, the scheduler is reset for each rollout step, reaching its starting learning rate after 10 warm-up iterations. We calculate the loss aggregated across the rollout steps, as in <cit.>.





 §.§ Hardware


The ddm was trained on compute nodes equipped with four AMD Instinct MI250X gpus, each with 128 GB memory. To speed up training, we used data parallelism across 32 model instances for stages A–C and 16 for stage D. For stages B–D, we used model parallelism to run one instance of the model on each node. This splits the nodes and edges of the graphs across multiple GPUs and is necessary for fitting the 31km global and 2.5km regional stretched-grid model within the available gpu memory. The model parallel approach is described in <cit.>.



§ EVALUATION



This section aims to evaluate the quality of the forecasts from the perspective of a public weather forecast provider. A majority of the forecast parameters we provide on our public weather forecasting app Yr, are deterministic. The key properties we seek for our forecasts are accuracy, by low squared error, and that the frequency distribution of the forecasts match that of the observations. The latter is important to ensure that aggregated statistic over longer time periods (e.g. maximum daily wind speed) is accurate, and are important parameters for users that are not sensitive to the exact timing of weather events. As these properties are often in conflict with each other, we evaluate both separately.

To evaluate the ddm from <ref> against competing nwp systems, we reserved a separate test period that we did not use when selecting the model configuration. The test period spans the time period from June 1, 2023 to May 31, 2024. The final model in Section <ref> was retrained on the combined training and validation period from Feb 6, 2020 to May 31, 2023 in order to increase the training period.



 §.§ Baseline NWP models


The ddm is evaluated against state-of-the-art nwp models used at MET Norway. Operationally, MET Norway relies on meps to provide high-resolution short-range weather forecasts for lead times up to 61 hours. This system is used for both automated weather forecasts for the general public and public weather warnings issued by duty forecasters. Additionally, they are used by downstream users such as energy, transportation, and agriculture. We extracted the control run and the ensemble mean from this modelling system.

We also include the control run from ifs at 0.1^∘ resolution, as this is a model we use to assess the added value of high-resolution models. 



 §.§ SYNOP observation network



The models are evaluated against observations from MET Norway's network of  synop stations. To allow for a fair comparison of models at different resolutions, the gridded data is interpolated bilinearly to the station locations. For temperature, both the temperature and the model terrain height are interpolated to the station point. The temperature is then adjusted based on a lapse rate of 6.5 ^∘C/km applied to the difference between the station elevation and the bilinearly interpolated model terrain height.



 §.§ Overall evaluation


The parameters we investigate are t, ws, and p6h. Additionally, we are interested in 24h aggregations, including daily minimum and maximum temperature, daily maximum wind speed, and 24h accumulated precipitation. Such aggregated values are used to summarize a day and are heavily used by our users. We also look at mslp as a diagnostic for the model's ability to capture general synoptic developments.





<ref> demonstrates how the stretched grid approach is able to seamlessly move weather systems from the global into the limited area domain. This case shows the storm Ingunn for a forecast initialized at 06Z on January 27th, 2024. Ingunn can clearly be seen to move into and develop within the limited area domain. Overall, the DDM appears in this case to have a similar capability to the current operational IFS and MEPS system in simulating large-scale features of wind speed and mslp. A notable difference is the slight spatial smoothing of the data-driven fields. 

<ref> shows a forecast from the same period zoomed in on the mountains in northern Norway. Like MEPS, the DDM also simulates strong winds on the mountain tops and weaker winds in the valleys. This is in contrast to the coarse resolution IFS, which fails entirely to represent the stronger winds in the mountains. This indicates that the DDM, despite some spatial smoothing, is able to represent many of the systematic wind speed features in mountainous areas.

In terms of rmse, the ddm outperforms both nwp systems for temperature (<ref>a) and 6h precipitation (<ref>c). For wind speed, the model performs similar to the MEPS control, but worse than the MEPS ensemble mean (<ref>b). This is due to the fact that the ddm has been trained on analyses from the control run, which tries to represent meteorological features that are unpredictable. The ensemble mean is significantly smoother and therefore scores better for all lead times. For mslp, the ddm is comparable to both nwp systems for lead times below 24 hours (<ref>d). However, the error growth is larger, and it is outperformed by both nwp systems for longer lead times. We found that mslp is sensitive to the area weighting of the regional domain, with a lower weight giving better mslp scores but worse for the other parameters. This could indicate some over-fitting, which we will investigate further. 





 §.§ Temperature


t is greatly affected by the complex topography and coastline of Scandinavia. This makes it a challenging parameter to predict.

For instantaneous temperatures, the ddm outperforms meps for RMSE (<ref>a). The ddm improves the rmse by around 24 hours compared to meps control and ensemble mean. Improvements of this magnitude usually require many years of model development. A major contribution to this improvement likely comes from the way meps assimilates surface temperature observations and carries that state over from the analysis into the forecast. This is evident in the large increase in error from lead time 0h to lead time 6h. meps assimilates the same observations we use to verify the forecast with. This leads to a low RMSE score for lead time 0h. However, a large part of this is lost already at lead time 6h. The ddm does not suffer from this problem because it is trained against analyses that have these observations assimilated and is better able to propagate information from the assimilated state forward in time.

As an example, <ref> shows that when the analysis increment is large, and the NWP model is not able to carry the increment forward, a large systematic error is present for the remainder of the forecast. Note that we verify the forecasts using the same observations employed in the data assimilation. This implies that we cannot expect to get the same forecast improvements over areas where no observations are assimilated. 

The ddm underestimates the standard deviation of temperatures across a 24-hour period (not shown). The consequence of this is that daily minimum temperatures are overestimated and daily maximum temperatures are underestimated. Despite these biases, the ddm outperforms meps in RMSE for both aggregated variables (<ref>a–b).







 §.§ Wind speed




To assess the model's ability to predict extreme events, we use the ets, also known as the Gilbert skill score. For a particular wind speed threshold, the frequency of forecasts and observations exceeding or not exceeding the threshold are computed, thereby classifying events into hits (a), false alarms (b), missed events (c), and correct rejections (d). The ets for a particular wind speed threshold is defined by:

    ETS = a-â/a+b+c-â,

where

    â = (a+b)(a + c)/a+b+c+d ;

The ets penalises forecasts with a high number of false alarms and event misses.

Spatial smoothing is significant for wind speed because of the double-penalty problem associated with mse loss. When a model predicts an event slightly off its actual position, it gets penalized twice: once for predicting the event at the incorrect location, and again for not predicting the event at the correct location. Smoothing the wind field improves the rmse but decreases the forecast's effectiveness in warning against strong wind events. The ddm outperforms the IFS in terms of ETS for most thresholds, which is largely due to IFS's poor ability to represent strong winds in the Norwegian mountains. It also has ETS scores close to meps, except for wind speeds above 15m/s (<ref>a). We also see that the ddm produces fewer strong wind events than what is observed (<ref>c). This is problematic for public weather forecasting. A simple solution is to adjust the distribution of winds to better match the distribution of climatology by a separate post-processing step.

For weather warnings, the timing of an event is often less important. To evaluate this need, we computed the 24h maximum wind speed (from the four instantaneous values within the time period). The meps control outperforms the ddm for an even greater range of wind speed thresholds (<ref>b) than for instantaneous wind speed. This is because the ddm has an increasing number of misses for higher thresholds, and has skill closer to the ensemble mean.





 §.§ Precipitation


For 6h accumulated precipitation, the ddm performs similar with respect to ETS to the NWP baseline models (<ref>a). As with wind speed, we see that the ddm underestimates higher precipitation amounts (<ref>c), which can be attributed to spatial smoothing.

Many of our users are interested in aggregated precipitation over longer periods, particularly at a daily timescale. Although the ddm is not directly trained to produce accurate 24h precipitation accumulations, we see that the model is competitive compared to the baseline models (<ref>b). The underestimation of extreme precipitation is less evident at the 24h timescale (<ref>d) than at the 6h time scale, and is likely due to the fact that any temporal smoothing at the 6h time scale does not negatively affect the 24h aggregated values.













For public weather forecasting, discriminating between a dry period and a rainy period is important as this affects the weather symbol we use to represent the period. For this, we use a threshold of 0.5mm, above which rain will appear in our weather symbol. <ref> shows that the ddm is better able to discriminate between the occurrence and non-occurrence of precipitation than the NWP models for lead times beyond 6h.



§ CONCLUSIONS AND FUTURE WORK



We have shown that a stretched-grid approach to regional data-driven modelling outperforms a state-of-the-art high-resolution NWP model for certain key parameters used in public weather forecasting, including instantaneous 2m temperature, 24h minimum and maximum temperature, and whether 6h precipitation exceeds 0.5mm. The model also provides competitive forecasts of instantaneous 10m wind speed, 24h maximum wind speed, and 6h and 24h accumulated precipitation, though the model tends to underestimate extremes.

The performance of the ddm is promising enough that it warrants use in public weather forecasting. However, several challenges must be overcome before we can use the model operationally. Firstly, forecasts at an hourly time scale is needed to serve the needs of our users. This will pose an even bigger challenge regarding GPU memory and it is still an open question if such a model can perform as well as a model trained with a 6h time step. Alternatives to auto-regressive approaches could be considered, such as training a separate model that interpolates 6-hour scenarios in time.

A second challenge is the need to provide probabilistic forecast <cit.> parameters, such as the quantiles of the probability distribution, which we use in our app to indicate the risk of heavy rainfall and strong winds. This can be solved by training an ensemble model or by directly modelling the quantiles using an appropriate loss function for scoring quantile forecasts.

To further improve the model in general, using more observational data could be important. For temperature, the ddm was able to utilize the information from assimilated observations better than the NWP model. Additional observations could be assimilated into the input grid before training, but a more general approach would be to use them directly with a separate encoder. Using crowdsourced observations, which have previously been shown to add value in operational weather prediction <cit.>, is an avenue we will explore further.

Due to their low computational requirements, ddms will allow us to run models for longer lead times than is typically affordable with high-resolution nwp models. This potentially allows us to use a single model for a range of timescales, instead of merging model runs from short-range, medium-range, and extended-range NWP systems as we do today. Operationally, these systems require separate post-processing as the models have different resolution and biases. End users would benefit from seamless forecasts without noticeable jumps across timescales.










 §.§ Acknowledgments

This work was supported by computing and storage resources provided by Sigma2 – the National Infrastructure for High-Performance Computing and Data Storage in Norway (Grant No. NN10090K).



 §.§ Data availability statement

The verification data (observations and forecasts) in this article are made accessible in Verif-format <cit.> on Zenodo (<https://zenodo.org/communities/verif/>).












ieeetr


  









