Real-time multi-model decadal climate predictions: Supplementary Information

Doug M. Smith1*, Adam A. Scaife1, George J. Boer2, Mihaela Caian3, Francisco J. Doblas-Reyes4, Virginie Guemas4, Ed Hawkins5, Wilco Hazeleger6,13, Leon Hermanson1, Chun Kit Ho5, Masayoshi Ishii7, Viatcheslav Kharin2, Masahide Kimoto8, Ben Kirtman9, Judith Lean10, Daniela Matei11, William J. Merryfield2, Wolfgang A. Müller11,Holger Pohlmann11, Anthony Rosati12, Bert Wouters6and Klaus Wyser3

1Met Office Hadley Centre, FitzRoy Road, Exeter, EX1 3PB, UK.

2Canadian Centre for Climate Modelling and Analysis, Environment Canada, Victoria, British Columbia, Canada

3Rossby Centre, Swedish Meteorological and Hydrological Institute, 60176 Norrköping, Sweden

4Institut Català de Ciències del Clima, Carrer del Doctor Trueta, 203 08005 Barcelona, Spain

5NCAS-Climate, Department of Meteorology, University of Reading, Reading, RG6 6BB. UK

6Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands

7Meteorological Research Institute, Japan Meteorological Agency, Tsukuba, 305-0052 Japan

8Atmosphere and Ocean Research Institute, University of Tokyo, Kashiwa, 277-8568 Japan

9RSMAS/MPO, University of Miami, 4600 Rickenbacker Causeway, Miami, FL 33149, USA

10Space Science Division, Naval Research Laboratory, Washington, D. C., 20375 USA

11Max-Planck-Institut für Meteorologie, Bundesstraße 53, 20146 Hamburg, Germany.

12Geophysical Fluid Dynamics Laboratory, Princeton University, Princeton, New Jersey, USA

13WageningenUniversity, Wageningen, The Netherlands

Key Words: Decadal climateprediction; Multi-Model ensemble; forecast

*Corresponding author: Doug Smith, Met Office Hadley Centre, FitzRoy Road, Exeter, EX1 3PB, UK,

Details of prediction systems

CCCMA

CCCma decadal forecasts are based on the CanCM4 climate model (Merryfield et al. 2011), which is similar to the CanESM2 model used for CMIP5 long-term projections except that the latter includes terrestrial and ocean ecosystem components. Atmospheric resolution is approximately 2.8 in latitude and longitude with 35 layers, whereas ocean resolution is approximately 0.94 in latitude by 1.4 in longitude with 40 levels. The forecasts are initialized from a set of assimilation runs, one for each ensemble member, begun in 1958 from different initial conditions drawn from a multicentury spinup run. Model variables in these runs are constrained toward full-field gridded observational atmospheric temperature, horizontal wind components and specific humidity (24 hour time scale), as well as sea surface temperature and sea ice concentration (3 day time scale). Data sources are ERA40 and ERA Interim (atmosphere), ERSST and OISST (sea surface temperature) and HadISST1.1 (sea ice), transitioning to Canadian Meteorological Centre equivalents in 2010. Subsurface ocean temperatures from the GODAS analysis are assimilated off-line prior to the forecast via a procedure similar to that of Tang et al. (2004), following which salinity is adjusted to prevent static instability as in Troccoli et al. (2002). Forecast biases in long-term temperature trends are removed as in Kharin et al. (2012) and aspects of historical skill are analyzed in Fyfe et al. (2012) and Boer et al. (2012).

GFDL

The GFDL assimilation system is based on the ensemble adjustment Kalman filter (EaKf; Anderson 2001), which is a deterministic variant of the ensemble Kalman filter. The EaKf estimates the probability distribution function (PDF) of climate states by combining the prior PDF derived from model dynamics and the observational PDF. It uses a two-step data assimilation procedure(the first step computes ensemble increments at an observation location and the second step distributes the increments over the impacted grids) for an ensemble Kalman filter under a local least squares framework. In both steps the filtering process is implemented by a multivariate linear regression with consideration of temperature and salinity covariance (Anderson 2003). The data adjusted ensemble members are the realizations of the analysis PDF and serve as the initial conditions for the next ensemble integration.

The GFDL system consists of an EaKf applied to GFDL's fully coupled climate model CM2.1 (Zhang et al. 2007), which is designed to produce a better balanced initialization as opposed to each component model using its own assimilation system. The ocean component of the ensemble coupled data assimilation (ECDA3.1) is the Modular Ocean Model version 4 (MOM4) configured with 50 vertical levels and 1° horizontal resolution, telescoping to 1/3° meridional spacing near the equator. The atmospheric component has a resolution of 2.5o x 2 o with 24 vertical levels. The first guess is given by a fully coupled model where the atmosphere is constrained by an existing atmospheric analysis. Ocean observations of temperature, salinity, and SST are assimilated using covariance structures from the coupled model. Argo observations are included as they became available in the post-2000 period. The cross-interface covariance structures in the GFDL system allow for fully coupled assimilation. For the ocean component, observed temperature and salinity profiles and SST are assimilated daily (see details in The atmosphere is constrained by an existing atmospheric analysis. All ECDA3.1 experiments are performed with a 12-member ensemble that is used to compute state estimation, ensemble mean, and the spread of the estimate. The ECDA3.1 also uses covariance inflation that is designed to enhance the consistency of upper and deep ocean adjustments, based on climatological standard deviation being updated by observations (Zhang and Rosati 2010).

Using the ECDA3.1 system, analysis covering the period 1960-2011 was obtained and from the ensemble members of the analysis coupled model initial conditions are produced. Predictions consisting of 10 ensemble members were run for 10 years beginning on 1 January for every year from 1960 through 2011, for a total of 5,100 years of integration. The CMIP5 historical radiative forcing used GHG, solar, volcano, aerosol for the 1960-2005 period and RCP4.5 scenario settings for predictions after 2005. A companion ‘uninitialized’ run, with the same forcings, was also made with 10 ensemble members from 1860-2040.

IC3/KNMI

The EC-Earth V2.3 has been used in this study. The EC-Earth V2.2 model and its main characteristics are described by Hazeleger et al (2010, 2012). In EC-Earth V2.3 a slightly different aerosol forcing has been used, consistent with the CMIP5 protocol. We use a horizontal spectral resolution of T159 (triangular truncation at wavenumber 159) and 62 layers in the vertical up to 5 hPa. The atmosphere model is derived from the Integrated Forecast System cycle 31r1 of the European Centre for Medium-Range Weather Forecasts (ECMWF). The ocean model is the NEMO version 2 model (Madec 2008) and the sea ice model is the LIM version 2 model (Goosse and Fichefet 1999). For details and further references we refer to Hazeleger et al (2012).

The system employs a full-field initialization. The ocean initial conditions have been produced with NEMOVAR at ECMWF, a multivariate 3D-var data assimilation method for the NEMO ocean model (Mogensen et al 2012). Observed 3-dimensional temperature and salinity and the sea surface height is assimilated. In particular, the NEMOVAR-ORAS4 five-member ensemble has been used, which is the operational analysis for the new Seasonal forecast system (S4) at ECMWF (Mogensen et al. 2012). The sea ice initial conditions come from a NEMO3.2-LIM2 simulation forced with ERAinterim (Dee et al, 2011) and nudged toward NEMOVAR-S4 for 5 members and from the GLORYS1V2 reanalysis for the 5 others. The atmosphere is initialized from ERA40 data (Uppala et al 2005) before 1989 and ERA-interim (Dee et al 2011) thereafter. The atmosphere is perturbed using singular vectors to create an ensemble of 5 members(see Du et al (2012) for further details).

MIROC5

Asystem for decadal prediction is based on MIROC version 5 (MIROC5) developed by University of Tokyo, National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology (Watanabe et al. 2010). The resolution of the atmospheric component is 1.41° of longitude and latitude with 40 vertical levels, and that for the oceans is 1.41°of longitude and 0.80° of latitude with 50 vertical levels. The model is assimilated from 1945 to 2011 with monthly gridded data of ocean subsurface temperature and salinity, interpolated daily. Here, temperature biases in historical expendable bathythermograph observations are removed (Ishii and Kimoto 2009). The model is forced by time-varying greenhouse gases and anthropogenic and natural aerosol concentrations. To avoid a model climate drift in decadal prediction, only observed anomalies are input for the initialization of MIROC5 (Tatebe et al. 2012). The model climatology is defined from an uninitialized model integration with the external forcing as a mean state from 1961 to 2000. The same period is used for the observational climatology. The assimilation scheme is based on an incremental analysis update, in which the model is relaxed to the observations being given increments proportional to differences between observed and model-predicted anomalies on a daily basis. A mixture of ensemble data assimilation and a lagged average forecasting is adopted for evaluation of the uncertainty of the decadal prediction. Three assimilation runs start in 1945 from different initial conditions taken from corresponding uninitialized runs integrated since 1850. There is no intercommunication among the three assimilation runs as an ensemble Kalman filter. A set of decadal predictions is conducted from initial states of the three assimilation models on July 1st and October 1st of the previous year and on January 1st of the first year of the individual prediction. Therefore the ensemble size is nine.

MOHC

The Met Office Hadley Centre Decadal Prediction System (DePreSys) (Smith et al. 2007, 2010) is based on the third Hadley Centre climate model (HadCM3) (Gordon et al. 2000), with a resolution of 2.5° x 3.75° in the atmosphere and 1.25° in the ocean. In order to create initial conditions for hindcasts and forecasts, HadCM3 is run in assimilation mode from December 1959 to the present day, including time-varying radiative forcing from changes in well-mixed trace gases, ozone, sulphate and volcanic aerosol, and solar irradiance. During this integration, the atmosphere winds, potential temperature and surface pressure, and the ocean temperature and salinity are relaxed towards atmospheric and ocean analyses, wherein values are assimilated as anomalies with respect to the model climate in order to minimize climate drift after the assimilation is switched off. The climatological period from which anomalies are computed is 1958 to 2001 for the atmosphere and 1951 to 2006 for the ocean. Atmospheric analyses are taken from ERA-40 and ECMWF operational analyses, while analyses of ocean anomalies are created using an updated version of the scheme developed by Smith and Murphy (2007), based on anomaly covariances calculated from HadCM3, with adjustments to improve the fit to observations. The relaxation timescales are 3 hours for the atmosphere and 6 hours for the ocean. Ensemble members are created by adding uncorrelated random perturbations (with standard deviation 0.005K) to sea surface temperature.

MPI

The MPI decadal climate predictions are performed with the system which was investigated for the decadal hindcasts undertaken for CMIP5 (Müller et al. 2012). For the decadal hindcasts and forecasts the MPI-M earth system model (MPI-ESM) is used which consists of the atmosphere component ECHAM6, land surface component JSBACH, ocean component MPI-OM and ocean biogeochemistry component HAMOCC. The components are coupled with the OASIS coupler. In a first step 3-dimensional oceanic temperature (T) and salinity (S) fields are produced by forcing the ocean component MPI-OM with surface fluxes from the NCEP/NCAR reanalysis. In a second step anomalies of these T and S fields are added to the MPI-ESM climatology and then assimilated (nudged) into the coupled model providing the initial conditions for the decadal predictions. In a third step these initial conditions are used to start the decadal predictions with the MPI-ESM (see Matei et al. 2012 for details). Ensembles of 10 members each are started on the 1stJanuary 2011 and 2012 using consecutive days as initial conditions.

MRI

The same methodology for decadal prediction as in the case of MIROC5 is used, except the climate model is MRI-CGCM version 3 developed by the Meteorological Research Institute (Yukimoto et al. 2012). The resolution of the atmospheric component is 1.125° of longitude and latitude with 48 vertical levels, and that for the oceans is 1.0° of longitude and 0.46° of latitude with 51 vertical levels. With this model, initialized and uninitialized model integrations are conducted and nine decadal prediction members are made following the CMIP5 protocol.

NRL

The Naval Research Laboratory has developed a linear climate model (NRL LCM) for estimating global and regional surface temperature from known sources of climate variability (Lean and Rind, 2008, 2009; Kopp and Lean, 2011). The model is constructed using linear multiple regression analysis of monthly averaged temperature observations (hadcrut3v) from 1978 to 2011 with ENSO (at three different lags), volcanic aerosols (at two different lags), total solar irradiance (at one lag), anthropogenic factors (RCP4.5), an annual cycle (2 terms), a semiannual cycle (2 terms), the QBO (at one lag) and NAO (at one lag). The model’s geographical grid is that of the surface temperature dataset (72 longitudes, 36 latitudes) and there is no interpolation to grid points that lack data. The small AO and SAO cycles are present presumably because of imperfect deseasonalization of the dataset. The NAO and QBO terms are included for completeness, but are sufficiently small that their omission has minimal effect on estimated global temperature. Separate versions of NRL LCM are also constructed for the GISS Land plus ocean and land only temperature datasets.

Forecasts of monthly averaged surface temperatures (globally and regionally) are made by inputting to the NRL LCM estimates of future solar irradiance (repeating cycle 23, with a period of 12 years), anthropogenic influence (RCP4.5) and annual and semiannual cycles, with the ENSO, volcanic aerosols, QBO and NAO influences remaining constant. Annually averaged surface temperature forecasts are determined as the average of the monthly forecast for each year, then adjusted to a baseline of 1971-2000 by subtracting the model’s average value for 1981-2000, adding the observed average value for 1981-2000 and subtracting the observed average value for 1971-2000.

The forecasts assume no net effect of either ENSO or volcanic aerosols. If, as is likely, this is not the case, then an adjustment can be made to the annually averaged surface temperatures. For example, ENSO activity is estimated to alter the forecast global surface temperature by TENSO=a  b  ENSO index (where the annual mean ENSO index varies from -1.3 to 1.58 and has a 1 range of 0.758). For the hadcrut3v dataset the scaling coefficients that convert the ENSO index to an equivalent global surface temperature anomaly (based on annual means from 1980 to 2011) are a = 0.247 and b = 0.080. The 1 estimate of the ENSO influence on annual global surface temperature is 0.3, which indicates that there is a 66% chance that the actual annual mean global surface temperature will be within 0.3 of the value forecast for a given year, providing there is no volcanic activity. Analogous adjustments can be made for volcanic activity.

RSMAS

The coupled general circulation model used for the RSMAS forecast system is the National Center for Atmospheric Research Community Climate System Model version 3 (CCSM3, Collins et al. 2006). The ocean initialization follows Kirtman and Min (2009) in that the Geophysical Fluid Dynamics Laboratory (GFDL) optimal interpolation ocean-analysis (Derber and Rosati 1989) was simply interpolated to the CCSM3 ocean model (POP) grid. There was no special treatment of the momentum terms. The three ensemble members are generated by using different atmospheric states, but the same ocean state (see Kirtman and Min 2009 for details).

SMHI

The SMHI prediction system uses the EC-Earth coupled model (Hazeleger et al, 2010). This is based on IFS atmospheric model (Bechtold at al, 2008) coupled to Opa-Nemo/Lim2 ocean/sea-ice model (Madec, 2008 ; Fichefet and Morales Maqueda, 1997). Atmospheric model resolution in these experiments is T159 for the spectral computations and an equivalent (1.125x1.125 deg) on Gaussian reduced grid for physics, with 62 sigma levels. The ocean and sea-ice models use a tri-polar grid (ORCA 1) with approximately 1x1 deg resolution and 42 z-levels.

The initialization uses the anomaly method for ocean and sea-ice thermodynamic variables. Ocean: 3 dimensional temperature, salinity and currents, sea-ice and snow-ice depths. Sea-ice: extent, surface temperature and velocity. We project the initial state onto the model's space by adding the initial daily anomalies (computed with respect to the observed climatology) to the coupled model climatology (both over 1960-2005). The ocean initial states and observed climatology are NEMOVAR-S4 analysis (Balmaseda et al, 2013), while for the ice these are provided from an ocean run forced by atmospheric analysis. A physical consistency check and adjustment procedure is applied to the initial anomalies.Atmospheric initial conditions are ECMWF operational analysis. The coupling frequency is 3h and the atmospheric model time step is 1h.

The external forcing in these experiments and their uninitialized counterparts was according CMIP5 with historical time-varying values of aerosols, mixed gases and ozone, sulphates and solar irradiance used used until 2005 and RCP4.5 thereafter.The ensemble members were obtained through combined perturbations of ocean (perturbed Nemovar analysis), sea-ice (observed and forced run simulated initial conditions) and atmosphere (time lag).

Reading (AR1 and CA)

The statistical prediction system (Ho et al. 2012) predicts the annual mean sea surface temperatures (SSTs) for the coming decade, based on the global gridded observation data set HadISST (with resolution of 1degree) from 1871 to 2011 (Rayner et al. 2003). We first decompose the SST anomaly time series for each grid box into two components: radiatively forced trend and internal variability, as these are modelled separately. The forced trend component is estimated by a linear regression of the SST time series on the 1-year-lagged time series of historical observed equivalent global mean CO2 concentration (Meinshausen et al. 2011). The residuals of the regression are taken as the internal variability. This internal variability component is modelled in two ways: (1) a first-order autoregressive (AR1) model (for all grid boxes globally) and (2) a constructed analogue (CA; van den Dool 2007) model (for the Atlantic only). For AR1, the model is fitted to the internal variability time series from up to 89 years before the start year of the forecast for each individual grid box. As for CA, a single model is fitted to spatial fields of SST internal variability, such that we develop a weighted linear combination of historical spatial fields which is closest to the spatial field at the start year. The future spatial fields are predicted by carrying forward the estimated weights. The SST anomaly prediction is the sum of the predicted forced trend based on the equivalent global mean CO2 concentration for the RCP4.5 scenario and the predicted internal variability by either AR1 or CA.