Ensemble prediction using a new dataset of ECMWF initial states-OpenEnsemble 1.0

Ensemble prediction is an indispensable tool of modern numerical weather prediction (NWP). Due to its complex data flow, global medium-range ensemble prediction has so far remained exclusively as a duty of operational weather agencies. It has been very hard for academia therefore to be able to contribute to this important branch of NWP research using realistic weather models. In order to open up the ensemble prediction research for a wider research community, we have recreated all 50+1 operational IFS ensemble initial states for OpenIFS CY43R3. The dataset (OpenEnsemble 1.0) is available for use 5 under a Creative Commons license and is downloadable from an https-server. The dataset covers one year (December 2016 to November 2017) twice daily. Downloads in three model resolutions (TL159, TL399 and TL639) are available to cover different research needs. An open-source workflow manager, called OpenEPS, is presented here and used to launch ensemble forecast experiments from the perturbed initial conditions. The deterministic and probabilistic forecast skill of OpenIFS (cycle 40R1) using this new set of initial states is comprehensively evaluated. In addition, we present a case study of typhoon Damrey from 10 year 2017 to illustrate the new potential of being able to run ensemble forecasts outside major global weather forecasting centres.

3 Initial state perturbations in the ECMWF ensemble

Singular Vectors
Singular Vectors represent the fastest-growing perturbations to a weather forecast -called the trajectory -within a finite time 85 interval (Lorenz, 1965;Buizza, 1994;Palmer et al., 1998). In order to compute singular vectors one linearises the governing equations around a given trajectory. The idea behind using singular vector based initial perturbations is that these are the dynamically most relevant structures and hence will dominate the forecast uncertainty (Ehrendorfer and Tribbia, 1997;Leutbecher and Palmer, 2008;Leutbecher and Lang, 2014). Here, growth is defined with respect to a specific metric. The metric used at ECMWF is the so called dry total energy norm (Leutbecher and Palmer, 2008) and singular vectors are computed with 90 an optimisation interval of 48 h and a spatial resolution of T L 42 and 91 vertical levels. Different sets of singular vectors are computed, the leading 50 singular vectors for the Northern and Southern Hemisphere, and the leading 5 singular vectors for each active tropical cyclone (see Buizza, 1994;Puri et al., 2001;Barkmeijer et al., 2001;Leutbecher and Palmer, 2008, for details). While singular vectors targeted at the extra-tropics are optimised for the whole troposphere, singular vectors targeted at tropical cyclones are optimised for growth below 500 hPa.

Ensemble of Data assimilations
At ECMWF an ensemble of 4D-Var data assimilations (EDA, Buizza et al., 2008b) is run to provide uncertainty estimates for both the ensemble forecasts and the high resolution analysis. The EDA consists of a number of perturbed members and one unperturbed control member. Originally it was run with 10 perturbed members but this changed in 2013 when the number of members was increased to 25. Recently the number of members was increased again to 50 . In this study 100 we make use of the EDA configuration with 25 perturbed members. For each EDA member the observations and boundary conditions are perturbed and the SPPT scheme is used to simulate the impact of model error.

Construction of initial perturbations
Perturbations based on the EDA short-range background forecasts are combined with singular vector based perturbations to build the initial conditions for the ECMWF ensemble forecasts. EDA perturbations are derived by subtracting the mean of the 105 EDA from each EDA member, resulting in a total of 25 perturbations. In addition to the EDA perturbations, 25 singular vector based perturbations are generated.
For each singular vector perturbation an individual linear combination from all singular vectors is constructed. This is done by scaling the singular vectors by random numbers drawn from a multi-variate Gaussian distribution (see Leutbecher and Palmer, 2008;Leutbecher and Lang, 2014).

110
The perturbations are centred on the high resolution analysis and have a plus-minus symmetry, i.e. the initial conditions of the first perturbed ensemble forecast member are derived by adding the first EDA perturbation and the first singular vector perturbation to the high resolution analysis. The initial conditions of the second perturbed ensemble forecast member are then The most relevant changes between this model version and the currently operational ECMWF version (CY46R1) that affect the ensemble initial states arise from changes to the model uncertainty representations (see Lock et al., 2019) and from the removal of the plus-minus symmetry. On top of this list, a number of changes have been introduced to the IFS model, the IFS data-assimilation process and amount of observations used in data-assimilation. These all will naturally also affect the ensemble initial states, due to direct or indirect contributions.

130
The data is tarred, gzipped and packed such that a single tarz-file contains all initial states (control state and perturbed initial states) for a given date and time (00 and 12UTC). The files are in the form "YYYYMMDDHH.tarz". The data is furthermore arranged into separate directories for the three available resolutions:  Table 1 lists the types of files a single downloaded tarz-file consists of. In each file there are 50 of the pan, psu, pua and pert files, one for each perturbed initial state. The first four files in Table 1 form an unperturbed atmospheric initial state that is used to initialize a control forecast. ICMGGoifsINIUA-files contain model level variables on a Gaussian reduced grid.
ICMSHoifsINIT-files contain model level variables in spherical harmonics representation. ICMGGoifsINIT-files contain a 140 number of surface level variables on a Gaussian reduced grid. ICMCLoifsINIT-files contain climatological surface fields for albedo at different radiation wave lengths, leaf area indexes, soil temperature and sea ice area fraction. For detailed description of the contents of each file we refer the reader to Appendix A. On top of these four files, OpenIFS requires static climatological files for radiation calculations and various namelists describing the hybrid sigma coordinates etc. In order to utilize the wave model, a separate set of wave model input files and namelists is required. These can be obtained from OpenIFS-support team on request. And finally, the pert_${MEMBER} files contain raw Singular Vector perturbations. These can be used to decompose the initial state perturbations into SV and EDA parts.
These files can be used to form four different kinds of initial states for experimentations, listed in Table 2. Control forecast and SV+EDA perturbed forecasts can be initialized directly from the provided files. For EDA only or SV only perturbations some file manipulation is required. In order to get initial states with EDA only perturbations, pert-file needs to be subtracted 155 from pua-file. An example on how to do this using Climate Data Operators (CDO; Schulzweida, 2019) software is given in Appendix B. The file manipulation is always subtraction: the +/-symmetry is build into the files, i.e. pert_001 and pert_002 will be identical fields with different signs. For SV only perturbations, the control state should be used, and the pert-files added to the ICMSHoifsINIT-files. Again, the same procedure with grid point conversion and +/-symmetry applies here. The workflow shown in Appendix B can be applied here with minor modifications.
160 5 Running ensembles from the dataset

Workflow manager -OpenEPS
In order to run large ensembles of a forecast model, a workflow manager is essential. For this purpose, we use here a simple yet efficient software called OpenEPS. We want to emphasize that the provided dataset of initial states is not tied to this, or any, software. OpenEPS is mostly written in bash but utilizes GNU make to handle parallel job executions in HPC, Linux cluster 165 or laptop environments. The software is freely available under an Apache 2.0 license (https://github.com/pirkkao/OpenEPS). Instructions on how to use the software, as well as a few example cases, are provided with the software download. We will nonetheless provide a concise description of the OpenEPS software here as the workflow would be similar with any workflow manager. The general workflow is handled as follows: 0. Setup a computing environment specific file containing various architecture settings 5. Link or generate initial states for the next date 6. Go back to 3 until all dates have been cycled through Step 0 only needs to be completed once for each unique computing environment. In the current HPC implementation, OpenEPS reserves all the wanted computing resources at the same time. For example, if the user wants to run a total of 5 185 ensemble members, each to be executed concurrently and each using 20 cores would mean submitting a batch job requesting a reservation of 100 cores for a timeslot of N minutes. If the user would want to run 50 ensemble members in total following the previous setup, OpenEPS would compute these in 10 consecutive batches within the same batch reservation. Instead of reserving 100 cores for N minutes, the cores would now be reserved for N*10 minutes. One could naturally increase the amount of computing resources as well in order to keep the execution time to a minimum, i.e. reserve 100*10 cores for N 190 minutes.
It is also possible to do online post-processing within the workflow, i.e. run scripts to manipulate each model forecast after they are finished. Note that due to the nature of the HPC implementation this means all the reserved resources might be sitting idle while this is happening. Usually the computing resources required for model forecast calculations are much larger than those used in post-processing the output, so caution is advised here. We recommend using this option only when online post-195 processing is required as part of the workflow, e.g. in algorithmic model tuning. Also, since manipulation of the initial state files is resource demanding, it it is highly advisable that modifications to the initial states (separation of SV and EDA parts for example) are done as a separate task before the actual model integrations. Workflow for this is also supported in OpenEPS.

Forecast model setup
Although the initial states have been generated using IFS version CY43R3, we use here the OpenIFS version matching IFS 200 CY40R1 as the forecast model. This is due to practical reasons: at the time of writing this paper, the matching cycle for OpenIFS (CY43R3v1) was still in preparation. We foresee that the forecast model difference will affect the testing somewhat, mainly due to differing analysis and forecast biases, i.e. the analysis bias is affected by the model cycle and hence a CY43R3 analysis and CY43R3 forecast model will have more similar biases than a CY43R3 analysis and a CY40R1 forecast model. This will result in potentially better scores when analysis and forecast were created by the same cycle. The forecasting skill in early forecast 205 lead times might also be somewhat degraded due to a stronger-than-usual spin-up effect. However, we still feel confident that the forecast skill evaluation is of value despite the forecast model differences. A number of physical parametrization and model dynamics changes happened between CY40R1 and CY43R3, we refer interested readers to the ECMWF OpenIFS and IFS webpages (ECMWF, 2020a, 2019a, b).

Experiment setup 210
We have run a number of experiments in order to assess how well the initial state perturbations fare with OpenIFS, these are listed in Table 3. The experiments cover the three provided horizontal resolutions: T L 159 (∼ 120 km), T L 399 (∼ 50 km) and T L 639 (∼ 32 km). Also, the different initial perturbation types (SV and EDA) were tested separately in order to illustrate the efficiencies of the perturbations in generating ensemble spread with the OpenIFS setup. All of the experiments were run without any model uncertainty representations.
215 Running large numbers of ensemble forecasts requires a substantial amount of computational resources.  demonstrated that the number of ensemble start dates is much more important than the size of the ensemble in extracting the mean probabilistic skill of the system. Thus, we keep the number of start dates high, but decrease the ensemble size for the on average over many cases and within sampling uncertainty caused by a finite number of cases and ensemble members (see Leutbecher and Palmer (2008), for an in-depth discussion). We use here operational ECMWF analyses from the forecast period as the truth. The model output is truncated to a 1 • /1 • regular grid before any other post-processing is done.  A mathematical prerequisite for calculating fair-CRPS is that the ensemble members need to be exchangeable, which is not 260 fulfilled in the ECMWF ensemble under this initial state construction style, where the SV and EDA perturbations are used with +/-symmetry (see . Note, this is no longer the case from 2019 onwards in ECMWF operational ensemble configuration. Therefore, the smaller sized ensembles here are either constructed from odd or even ensemble members, which fulfills the prerequisite. As per construction, the normal CRPS scores the better the more ensemble members the ensemble contains (Fig 3). But, the fair-version of the score gives near identical results for the different ensemble sizes. This allows us 265 to meaningfully compare different configurations.    The SV perturbation amplitude is a tuning parameter of the ensemble (Leutbecher and Palmer, 2008;Leutbecher and Lang, 2014). Figure 5 illustrates the sensitivity of T L 159 resolution ensemble to a change in the initial perturbation amplitude. There  7 Applications in case studies: an example from forecasting a tropical cyclone The following case study is aimed to illustrate the usefulness of the new data set and the OpenEPS repository. We want to emphasize that the goal here is not to dissect how the ensemble behaved or what caused the differences in the forecasts, but to give an idea how running ensembles can provide plenty of insights normally unavailable from a single model forecast.
Typhoon Damrey started as a tropical depression in the Philippine archipelago on October 31st, 2017. After moving across 285 the open sea to the west of Philippines, it started to rapidly intensify and reached its peak strength on the third of November (the control forecast initial state for MSLP and 200 hPa wind vectors for November 2nd 12UTC is illustrated in Appendix C Fig   C1). The typhoon made a landfall in Vietnam the following day and caused severe damage and loss of life (see e.g. GFDRR, 2018). Report of the operational forecasting performance of the event can be found from the ECMWF Severe Event Catalogue (ECMWF, 2020b). 290 We use our dataset to launch a 20-member OpenIFS ensemble starting from November 2nd 12UTC with both SV and EDA perturbations active with T L 639 resolution. The ensemble mean MSLP and ensemble spread for the 0th time step of the ensemble forecast is shown in Fig 6. The observed track of Damrey is also plotted (red). Notably, the largest differences in the initial states are focused on the South side of the typhoon core (MSLP minimum). Large differences can also be observed in the East-West structure of the typhoon. Lang et al. (2012) show how especially the SV perturbations can rapidly alter the 295 TC location and intensity. Model forecast differences (due to the initial state perturbations) after a 36 h forecast lead time are shown in Fig 7. The ensemble mean as well as majority of the ensemble members place the landfall location too South and propagate the typhoon core too slowly (too Easterly location). The exact location and timing of the typhoon landfall is however within the likely solutions from the ensemble.
These kind of case studies using ensembles could be used to study various mechanics of the model: Was there already 300 something incorrect in the TC structure in the unperturbed initial state? Did some initial state perturbations correct this and made an impact on the TC forecast? Were local perturbations to the typhoon core essential, or was it perhaps the mean flow forcing that was more important to get correct? Would activation of the OpenIFS wave model improve all of the forecasts due to improved representation of momentum exchange between ocean surface and atmosphere?   For all resolutions, SV perturbations generate least spread in the ensemble. Nonetheless, having SV perturbations active on top of EDA perturbations clearly brings value to the forecast skill of the system. All perturbation types increase the accuracy of the ensemble mean when compared against a control forecast initialized from an unperturbed analysis state. We have also tested the impact of inflating the amplitudes of the initial state perturbations. Increasing the amplitudes of the initial state perturbations result in an increase of the forecast skill of the system.

320
Inspection of especially the lowest resolution experiments reveals that all the ensembles are underdispersive, i.e. the ensemble spread is much smaller than the ensemble mean RMSE. Also, the EDA perturbations in the Tropics are inefficient in generating spread in forecast lead times less then 96 h. This is expected since our ensemble configuration is missing a model