Clim Dyn (2009) 33:233–243 DOI 10.1007/s00382-008-0456-9
MVL spatiotemporal analysis for model intercomparison in EPS: application to the DEMETER multi-model ensemble J. Ferna´ndez Æ C. Primo Æ A. S. Cofin˜o Æ J. M. Gutie´rrez Æ M. A. Rodrı´guez
Received: 11 February 2008 / Accepted: 12 August 2008 / Published online: 27 August 2008 Ó Springer-Verlag 2008
Abstract In a recent paper, Gutie´rrez et al. (Nonlinear Process Geophys 15(1):109–114, 2008) introduced a new characterization of spatiotemporal error growth—the so called mean–variance logarithmic (MVL) diagram—and applied it to study ensemble prediction systems (EPS); in particular, they analyzed single-model ensembles obtained by perturbing the initial conditions. In the present work, the MVL diagram is applied to multi-model ensembles analyzing also the effect of model formulation differences. To this aim, the MVL diagram is systematically applied to the multi-model ensemble produced in the EU-funded DEMETER project. It is shown that the shared building blocks (atmospheric and ocean components) impose similar dynamics among different models and, thus, contribute to poorly sampling the model formulation uncertainty. This dynamical similarity should be taken into account, at least as a pre-screening process, before applying any objective weighting method. Keywords Ensemble prediction systems Error growth Model intercomparison MVL diagram Spatiotemporal dynamics
J. Ferna´ndez (&) A. S. Cofin˜o Department of Applied Mathematics and Computing Sciences, University of Cantabria, 39005 Santander, Spain e-mail:
[email protected] C. Primo European Centre for Medium-range Weather Forecasts, Reading, UK J. M. Gutie´rrez M. A. Rodrı´guez Instituto de Fı´sica de Cantabria, CSIC-UC, Santander, Spain
1 Introduction Ensemble prediction systems (EPS) are the most widely used methods to assess the uncertainties inherent to the meteorological prediction at different time scales (see, e.g., Gneiting and Raftery 2005; Palmer 2002). EPS exploit the use of multiple realizations (an ensemble) of a prediction system with varying initial conditions, parameters and/or models, depending on the uncertainties to be evaluated. The results are analyzed working in a probabilistic framework, usually assuming a prior uniform distribution for the realizations (equal weights). In the case of multimodel ensembles, alternative weighting schemes have been proposed in the literature, such as the reliability ensemble averaging (REA method, see Giorgi and Mearns 2003), the Bayesian model averaging (BMA, Raftery et al. 2005) or the Forecast assimilation (Stephenson et al. 2005), since the assumption of equiprobable models is not justified in practice; these weighting methods are based on the expected performance of model output variables (model skill) and, thus, those models better representing the climate are favored when combining the results. However, this approach provides only a partial solution to the weighting problem, since it neglects the similarity of the existing models, which are based on a reduced number of building blocks developed at major research centers (ocean model, atmosphere model, couplers, etc.), with preferred combinations. Therefore, in ensembles with favored models or building blocks, a first intercomparison analysis regarding model similarity would be necessary in order to obtain proper results. In this paper we analyze the problem of model intercomparison considering a recently introduced technique: the mean–variance logarithmic (MVL) diagram (Gutie´rrez et al. 2008). This diagram has a solid theoretical basis
123
234
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
(Lo´pez et al. 2004; Primo et al. 2007) and allows characterizing EPS by analyzing the growth of the corresponding perturbations (e.g., introduced in the initial conditions for a particular spatiotemporal model). The MVL diagram is obtained representing the temporal and spatial components of the growth dynamics in the two axis of the MVL diagram, hence providing new insight into the interplay between these components, which is usually neglected in standard studies of spatiotemporal chaos. Gutie´rrez et al. (2008) applied this methodology to toy models (the Lorenz 96 system) and to Global Circulation Models (the monthly and seasonal ECMWF systems). They showed that, when applied to ensembles generated by perturbing the initial conditions, the MVL diagram exhibited a characteristic ‘‘fingerprint’’ governed by the ensemble generation procedure (random, bredding, etc.), in a first perturbationgrowing stage, and by the intrinsic model fluctuations, in a second stage after saturation. However, the behavior of the MVL diagram when applied to different GCMs was not addressed in the above study. The present work is a follow up of Gutie´rrez et al. (2008), where we extensively apply their methodology to the characterization of multi-model EPS, focusing on the intercomparison of the different models forming the ensemble from a dynamical point of view (dynamical similarity). This kind of multi-model ensembles have been studied, among others, in a series of EU funded seasonal forecast projects: PROVOST (Palmer and Shukla 2000), DEMETER (Palmer et al. 2004) and ENSEMBLES (Hewitt and Griggs 2004). These projects focus on the uncertainties associated with the initial conditions and the GCM formulation. In particular, DEMETER provides the longest seasonal ensemble prediction (seven models run in hindcast mode for 22 years) available to date. In this paper we consider this data set to apply the MVL approach. This paper is organized as follows. Section 2 describes the DEMETER simulations and analyzes the different ensemble generation procedures used in the project. Section 3 briefly introduces the MVL representation, which is later systematically applied to the DEMETER simulations in Sect. 4 obtaining the components with larger MVL variability. Section 5 describes the characterization of initial and stationary regimes in the MVL, associated with the ensemble generation procedure and the model dynamics, respectively. This section also presents some suggestions for model intercomparison using this approach. Some final remarks and conclusions are summarized in Sect. 6.
2 The DEMETER multi-model ensemble The EU-funded DEMETER project (Palmer et al. 2004) aimed at advancing in the concept of multi-model
123
ensemble prediction by performing seasonal simulations with seven different coupled atmosphere–ocean GCMs using nine different initial conditions. The experiments were performed in hindcast mode by initializing the models four times per year on 1st February, May, August and November and running for a period of 6 months. The seven models were run for different periods, but they share a common 22-year period from 1980 until 2001, which was used in this work. Table 1 shows the seven centers contributing simulations to the project from GCMs formed by different combinations of the atmospheric and oceanic components (the building blocks) shown in Table 2; note that four out of seven models include similar versions of the same ocean component (the OPA ocean model) and, thus, the seven models cannot be considered equiprobable in the resulting ensemble. In this paper we show that the MVL diagram can be considered a fingerprint of model dynamics and, therefore, it can help us to identify model similarities for pre-screening and quantification purposes. A common procedure to generate the perturbations of the sea surface temperature (SST) and wind stress was planned at the beginning of the project. According to Palmer et al. (2004), all models except the Max Plank Institute for Meteorology (smpi) used nine different initial conditions for the ocean prepared in a two-step approach. The starting point was a control ocean analysis forced with momentum, heat and mass flux data from the ERA40 reanalysis (Uppala et al. 2005); then, two perturbed ocean analysis were created by adding daily wind stress perturbations to the ERA40 momentum fluxes. In a second step, SST perturbations were added and subtracted from these ocean analyses in order to represent the uncertainty in SSTs, leading to a total of nine different initial conditions (considering also the unperturbed one). Atmospheric and land-surface initial conditions were taken directly from ERA-40. Thus, a total of 8 perturbed plus a control initial conditions were available at the start of the hindcast periods.
Table 1 Centers producing the GCM simulations for DEMETER and the corresponding atmospheric and oceanic components as labeled in Table 2 Label
Atm
Ocn
Simulation center
cnrm
A
a00
Me´te´oFrance
crfc
A
a
CERFACS
lody
B
a
LODYC
scnr
D0
a0
INGV
scwf
B
b
ECMWF
smpi ukmo
D C
d c
Max Plank Institute UK-MetOffice
The label field is used in all figures to identify the models
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS Table 2 Atmosphere and ocean model versions used in DEMETER, along with their horizontal resolution and vertical levels. Each model is labeled with a letter Resolution Atmosphere models
235
2.3 INGV The Italian National Institute of Geophysics and Vulcanology (scnr) produced the ocean initial conditions by the standard method after a 6-year spin up of the model with the ERA-15 fluxes. Unlike the standard procedure, atmospheric initial conditions were obtained from an AMIP-type run performed with the atmosphere component of the model forced with observed SSTs for the period 1972–2001.
A
ARPEGE
T63 L31
B C
IFS HadAM3
T95 L40 2.5° 9 3.75° L19
D
ECHAM5
T42 L19
D0
ECHAM4
T42 L19
2.4 UK-MetOffice
a
OPA8.2
2.0° 9 2.0° L31
a0
OPA8.1
2.0° 9 0.5°–1.5° L31
The UK-MetOffice (ukmo) followed the standard procedure described by Palmer et al. (2004).
Ocean models
a00
OPA8.0
182GP 9 152GP L31
b
HOPE-E
1.4° 9 0.3°–1.4° L29
c
GloSea
1.25° 9 0.3°–1.25° L40
d
MPIOM1
2.5° 9 0.5°–2.5° L23
Different versions of the same model are labeled with the same letter and a prime
However, each of the centers used slightly different perturbation methods; the final perturbations applied to the different models are not always clear from the available documentation.1 A brief description of the available information is given below. 2.1 ECMWF The ECMWF (scwf) followed the Palmer et al. (2004) procedure to generate the ocean analyses with their ocean model and added a strong relaxation to observed SSTs (Reynolds 2D-Var SST). They also changed the salinity beneath the surface layer to preserve the T–S relationship. Additionally, they added a weak relaxation to subsurface climatological salinity and temperature data and a geostrophic correction to velocity following a density change. The ocean model was initialized from an ocean at rest with January climatological temperature and salinity, and spun up for 4 years forced with monthly mean climatological wind stress, before starting the forcing with ERA40. 2.2 Me´te´o-France Me´te´o-France (cnrm) has followed the standard procedure with their ocean model (OPA), but starting three forced integrations at the beginning of each ERA40 stream (January 1958, January 1973, October 1986) from the same climatological oceanic state. This climatological state was obtained from an 8-year spin-up starting at rest and using ERA15 climatological monthly mean surface fluxes. 1
http://www.ecmwf.int/research/demeter.
2.5 CERFACS The European Center for Research and Advanced Training in Scientific Computation (crfc) spun up the ocean model from climatological values for temperature and salinity and from rest for the velocities, using a blended climatology of observed winds, ERA15 heat fluxes and gridded observed precipitation and evaporation data. After this 2-year spin up, during which the state remained close to the climatology, the model was forced with ERA40 daily fluxes, winds, and with a restoring SST term from November 1986 onwards. In order to produce perturbed ensembles of ocean initial conditions, an initial condition from the unperturbed forced ocean run was used every 3 months to start two wind-perturbed forced runs. Two weeks before the target date, four SST perturbations were linearly added and subtracted during 7 days to the restored SST and then persisted. This procedure was chosen in order to preserve the dynamical balance in the mixed layer scheme of the ocean model. The ocean initial conditions produced at CERFACS were shared with LODYC (see below). 2.6 Max-Plank Institute for Meteorology The MPI (smpi) hindcasts are initialized in a completely different way than the rest of the models. They used a coupled assimilation scheme in which the model is run in coupled mode up to the start of the hindcast with SST strongly relaxed to observations. They argue that this method provides a balanced set of atmospheric and oceanic initial conditions, but then the individual ensemble members are generated by choosing different atmospheric conditions around the hindcast start date. 2.7 LODYC The Laboratory for Dynamic Oceanography and Climatology (lody), used the ocean initial conditions produced by
123
236
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
CERFACS for the project since both share the ocean model. In summary, every coupled model used its own oceanic component to generate a set of 9 ocean initial conditions following similar practices. Then, the atmospheric and soil conditions from ERA-40 were used as atmospheric initial conditions for every member and model. The only exceptions to this procedure are the scnr model, that used an atmospheric state generated by its own atmospheric component forced with observed SST, and the smpi, that used a completely different method. With the exception of smpi, the atmospheric state at time 0 for all members of each model is identical. The perturbation starts from the surface over the oceans and we consider it from forecast day 2–90. Notice that this paper focuses on the characterization of the intrinsic spatiotemporal dynamical behavior of different GCMs and, thus, it is not concerned at all with how well each model represents the real climate. The assessment of how a perturbation is grown in the real climate system is not trivial and will be dealt with in a forthcoming paper. The reader interested in the performance of each model is referred to the Website of the project, where a comprehensive set of validation charts are presented.
3 The MVL representation The MVL diagram was introduced by Gutie´rrez et al. (2008) to characterize ensembles, based on the theory of spatiotemporal growth of perturbations in non-linear systems, adapted from the kinetic rough interface growth theory (see Lo´pez et al. 2004; Primo et al. 2007, 2008). This section summarizes, from a practical point of view, the computation and interpretation of the MVL diagram for the case of an ensemble with multiple models and multiple initial conditions. We also provide a brief physical interpretation of the two main quantities involved in the analysis, referred to as M value and V value (see below). For further theoretical background, the reader is referred to the above mentioned references. The MVL representation is based on the log-perturbations, defined as the logarithm of the absolute value of the differences between each of the single model ensemble members and their control run. The spatiotemporal evolution of perturbations for a single model is characterized by a curve in the MVL diagram with an initial stage driven by the growth of the initial perturbations and a final stage corresponding to the ‘‘climatological’’ model fluctuations, attained when perturbations saturate and can no longer grow. Therefore, this diagram can be considered a fingerprint of the ensemble and model dynamics. The selection of the control run is not relevant for the results so the
123
member labeled by each center as member 0 was used as control. For each of the m = 1,...,7 DEMETER models we define j = 1,...,8 log-perturbations for the members 1 through 8 as: m m hm j ðt; xÞ ¼ logj/j ðt; xÞ /0 ðt; xÞj;
ð1Þ
where /m j represents an undetermined variable of the member j of model m and /m 0 is, thus, the value of the same variable in the control run. The MVL diagram is a two-dimensional representation displaying in the abscissa the spatial average of the logperturbations (the M value): 1 X m Mjm ðtÞ ¼ h ðt; xÞ ð2Þ Nx x j and in the ordinate the spatial variance (the V value): 1 X m Vjm ðtÞ ¼ ðh ðt; xÞ Mjm ðtÞÞ2 : ð3Þ Nx x j The evolution of a single member simulation in the diagram is given by the parametric curve (M(t), V(t)). This two statistics were area weighted to account for the smaller spatial representation of grid-points near the poles (this has not much impact on the results, though). In the following we consider the averaged curves for the different members in order to obtain more representative results. It can be shown that the horizontal axis (M value) represents the spatially averaged error growth in time (given by the leading Lyapunov exponent), whereas the vertical evolution (V value) represents the spatial dynamics given by a changing spatial correlation length (characteristic spatial scale) of the perturbation field; the resulting spatiotemporal growth is a combined effect of both interplaying components (Gutie´rrez et al. 2008). Thus, a single-model EPS with perturbed initial conditions can be simply characterized using this diagram, which describes the corresponding growth of initial perturbations. On one hand, EPS systems evolve in time from left to right in the MVL diagram due to the exponential growth of perturbations; moreover, since non-linear systems are bounded, the perturbations saturate when their size reach the amplitude of the system, so the M value grows until saturation due to the ‘‘nonlinear barrier’’ (see Fig. 1). On the other hand, EPS systems evolve upwards (downwards) in the MVL diagram indicating a gain (loss) of spatial structure given by a growth (decay) of the spatial correlation. The vertical dynamics also saturates due to the finite spatial size of the system, since the correlation length cannot grow larger than the finite length of the system which imposes a ‘‘finite size barrier’’ (see Fig. 1). Thus, the V value dynamics is richer and more informative since it indicates the structural changes in the spatial properties of the growing
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
237
attractor), then, the spatial structure cannot grow larger and the V value remains stable until the non-linear barrier is reached and the MVL trajectory decays to the stationary point (line C). This kind of initial perturbed conditions are generated by special techniques such as Bred Vectors (Toth and Kalnay 1993) or Singular Vectors (Molteni et al. 1996) designed to keep the dynamical consistency.
4 DEMETER characterization
Fig. 1 Typical MVL trajectories for A random and spatially uncorrelated, B spatially correlated and non-assimilated, and C assimilated perturbations. See text for details
perturbations undergoing different regimes. The resulting saturation point (or area, for non-stationary models) corresponds to model ‘‘climatological fluctuations’’ (e.g., differences between two arbitrary states of the system), whereas the transient evolution is characteristic of both the ensemble generation procedure (assimilated and not assimilated initial conditions follow a different evolution in the diagram) and the model dynamics (location of the finite-size and non-linear barriers). For instance, Fig. 1 shows three typical trajectories in the MVL diagram for different initial perturbations. If the initial perturbation is random and spatially uncorrelated (line A), then the spatial correlation is small (for an uncorrelated and normally distributed initial perturbation V(0) = 1.23) and grows while it is assimilated into the system, gaining the system characteristic spatial structure until decay due to the non-linear effects. The initial V value is higher when the model is initially perturbed with a spatially organized pattern (Fig. 1, lines B and C). In such a case, there are two typical behaviors. On one hand, if the perturbed state is not assimilated into the current flow state (i.e., it is not in the attractor of the system), then the spatial structure of the initial perturbation is progressively destroyed in the dynamical assimilation procedure2 and, thus, the V value is first reduced (line B, in Fig. 1) and then increased when a new spatial structure, dynamically assimilated by the system, is created. The MVL trajectory then arrives at the same saturated point corresponding to climatological fluctuations, after reaching the finite size and non-linear barriers. On the other hand, if the initial spatial structure is flow-assimilated (i.e., it is on the system 2
Note the difference between the artificial assimilation procedure used in data assimilation (the result is something lying between a model state and the observations) and the dynamical assimilation of an initial condition lying out of the current flow state, which is dumped into the model attractor by the system dynamics.
The DEMETER multi-model seasonal hindcast offers a great opportunity for the comparison of the dynamics of the spatiotemporal growth of perturbations in different models. The high amount of degrees of freedom (7 models, 9 members, 4 seasons, 22 years, many variables, ...) limits the number of dimensions which can be simultaneously studied. Thus, in this section we analyze the variability of the different components in the MVL diagram and obtain the particular averaging displaying the maximum information and better characterizing model differences in the EPS dynamics. In particular, we restrict ourselves to the analysis of two variables: one close to the perturbed fields, surface temperature (T2m), and the other one far away from the perturbations, 500 mb geopotential height (z500). 4.1 Interannual variability of the MVL trajectories In the case of the DEMETER seasonal hindcasts, the time dependence spans three dimensions: Mjm ðtÞ ¼ Mjm ðf ; y; sÞ Vjm ðtÞ ¼ Vjm ðf ; y; sÞ;
ð4Þ
where the true time dimension is the forecast time (f) in days, while the year (y) and starting month or season (s) add dimensions which will allow us to analyse the interannual variability (in this section) and the seasonal dependence of the trajectories (in Sect. 4.2) on the MVL diagram. Figure 2 shows ensemble member averages for 2-m temperature for every model and for every year for the November initialization: D E D E Vjm ðf ; y; s ¼ NovÞ vs. Mjm ðf ; y; s ¼ NovÞ ð5Þ j
j
It is clearly apparent from this figure that the MVL representation is very sensitive to the model dynamics: Year to year differences are much smaller than model to model differences. This is a unique feature which cannot be achieved using directly the model output and focusing in the skill to reproduce a given variable. In a given year, all model simulations start from very similar initial conditions, and all models solve similar equations, so all models exhibit similar evolutions for a given variable; however these evolutions differ from year to year (which start from different initial
123
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
238
Fig. 2 MVL diagram of the 2-m temperature ensemble evolution for the November initializations. 22 curves (1 per year) plus the mean (in bold) are shown for each model and represent an average of the eight single-model ensemble perturbations. The dots on the lines represent forecast days. The first nine forecast days are labeled with their corresponding number
conditions); then, the interannual variability of a given variable is larger than the model-to-model variability. However, when we work with variable perturbations (differences between members) in the logarithmic space using the MVL diagram, we obtain information about the spatiotemporal dynamics of the model, which grows structures from the perturbed initial condition until saturation in the climatological fluctuation area. Thus, as it is shown in Fig. 2, interannual variations in the MVL diagram become negligible compared to model to model differences. Therefore, this diagram gives a fingerprint of the model dynamics which may provide a form of model intercomparison. It is, thus, reasonable to work with yearly averages, focusing on the model to model differences: D E V m ðf ; sÞ ¼ Vjm ðf ; y; s ¼ NovÞ jy D E ð6Þ Mm ðf ; sÞ ¼ Mjm ðf ; y; s ¼ NovÞ
Fig. 3 MVL diagram of the mean z500 ensemble evolution for the November initializations. Each curve is an average of the eight singlemodel ensemble perturbations and the 22 initializations for each model. The saturation area (within the outlined corners) is zoomed in Fig. 7
simulations starting in February, May and August. Figure 4 shows the seasonal dependence of the MVL diagram for a single model. The rest of the models exhibit a similar behavior. Again, this representation is mainly insensitive to the season and shows a similar evolution as that already shown for the November initialization. There appear only small differences in the mean V value and the final saturation region according to the seasonal cycle. However, model-to-model differences are still much larger than season-to-season differences. Thus, the rest of the study will focus on the November initialization, the conclusions being also valid for the other seasons.
jy
as illustrated also in Fig. 2 (bold lines with dots). For instance, lody and scwf nearly overlap in the diagram, indicating a similar dynamical behavior. After Table 1, we see that lody and scwf share the atmospheric component (IFS model). This similarity is also found analyzing the z500 variable (see Fig. 3) although, in this case, we observe a different initial transient characteristic of this variable (see Sect. 5.1). 4.2 Seasonal variability The previous results were shown only for the DEMETER simulations starting in November. There are three other
123
Fig. 4 MVL diagram of the mean ensemble evolution for the LODYC model. Each curve is an average of the 22 initializations for each starting month
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
239
Fig. 5 MVL diagram of the mean ensemble evolution for the c November initializations for 2 m temperature. Each model is shown separately and also the averages of the 8 single-model ensemble perturbations (left) and the 22 yearly initializations (right). The background ‘‘shadow’’ represents each single evolution (one member for 1 year)
4.3 Member variability The thinnest lines in Fig. 2 still show an average over the ensemble members (j = 1, ..., 8) for the 2-m temperature variable. Figure 5 decomposes, for each model, these lines into their eight members to assess the spread lost by ensemble averaging. The spread is still smaller than the intermodel spread and the ‘‘spaghettis’’ corresponding to each model evolve in a characteristic manner in the diagram (Fig. 6 shows the same analysis for z500). The yearly averages of single-model ensemble members (Figs. 5, 6, right) show, in general, less spread than the member average for each year (left panels) and provide information about the intra-model perturbation procedure and the model dynamics. The dispersion of the members is larger for the 2-m temperature, which is strongly affected by the perturbation procedure used to obtain the ensemble of initial conditions (this fact will be discussed later in the paper).
5 The MVL model fingerprint The different averaged views of the MVL trajectories introduced in the previous section allow us to examine the evolution of the different ensemble generation procedures (in the initial transient regime) and models (in the stationary regime). In this section we examine both regimes, which provide a fingerprint of the perturbation schemes and model dynamics, respectively. 5.1 Initial transient: perturbation schemes The first stage of an EPS trajectory in the MVL diagram, before reaching the nonlinear barrier, is determined by the particular perturbation scheme applied to the initial condition (see Fig. 1). In the case of the DEMETER multimodel ensemble, the initial perturbation was applied only on the surface, and the atmospheric variables considered were identical at time 0 (except for smpi). Therefore, the first point shown on all MVL diagrams correspond to day 2 (the MVL state is undefined if the difference between perturbed member and control is zero). We can clearly distinguish three different behaviors for the 2-m temperature case (see Fig. 2). First, most of the models start with a spatially coherent perturbation (large V value). However,
123
240
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
the structure is not compatible with the current state of the system and the V value decays until a new structure is created in the dynamical assimilation process. Second, there are two models, lody and scwf, that start from random spatially uncorrelated perturbations. The mean size of the log-perturbation (M value) at the start is an order of magnitude greater than in the previous models. Finally, the smpi model starts with saturated perturbations (this is clearer in Fig. 5). In this case, the model initial perturbations exhibit a size and spatial structure similar to the differences between two random states of the model (climatological fluctuations). Instead of the standard perturbation procedure, an ensemble of different initial conditions several days apart around the target one (lagged method) was used in this case. Two processes that take place simultaneously along the first simulated day, and before the first MVL point on the diagram can be computed, might be involved in the differences observed among models. On one hand, the SST and the atmospheric conditions were generated independently and are not coupled to each other. Thus, the atmospheric state needs to be adapted to the SST below. On the other hand, the atmospheric state comes from ERA40 and is not a proper state of the atmospheric model. Thus, this initial state needs to be driven towards the attractor of the corresponding model. The two models (scwf and lody) using the same model version (IFS, cycle 23r4) as that used to generate the ERA40 initial atmospheric conditions exhibit after the first simulation day near random and larger perturbations than the rest of the models and evolve according to typical MVL trajectory (Fig. 1a). For the other models (except smpi), the ERA-40 initial atmospheric condition induces a suddenly gained spatial structure which is not on the attractor and is destroyed before the perturbations can grow (following another usual behavior shown in Fig. 1b). Additional data of the first simulation day would be crucial to understand how this initial non-assimilated structure is gained. The results for the geopotential height (Fig. 3) are less responsive to the ensemble generation procedure, since the atmospheric state is the same for all simulations and the perturbations are introduced at the surface level (in the atmosphere–ocean coupling) and are seen from the 500 mb level as random perturbations. The similarities in the atmospheric component of the model mark the dynamical behavior of this variable. Those GCMs sharing the atmospheric component completely overlap their MVL trajectories. This is the case for the pairs lody–scwf and cnrm–crfc. The scnr and smpi models also share a very similar atmospheric component (based on different versions of the ECHAM model), however the ensemble generation procedure of smpi completely determines the MVL trajectory.
123
Fig. 6 As Fig. 5, but for 500 mb geopotential
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
When looking at the details of the different members of the ensemble (averaging the 22 years), shown in Figs. 5 and 6, we see that all members follow a similar MVL trajectory in the z500 case; however, in the case of 2-m temperature, some models exhibit a systematically different behavior for different members of the ensemble. For instance, the models from cnrm and ukmo show a different initial dynamics for two of the member initializations; however, after the initial transient all members follow a similar dynamics (that is the stationary regime characteristic of the model dynamics, which is similar in all members). This behavior could be related to the ensemble generation process. The nine members (including the control) are grouped in groups of three members; either sharing the same SST perturbation before the hindcast or the same wind stress perturbations when generating the ocean analyses. Thus, whichever member is taken as control member, there are always two other members sharing the same SST or wind stress perturbations and, therefore, two members differ only a little and the other six are further removed. On the other hand, crfc exhibits an extreme behavior where all the perturbed members seem to have a different dynamics all throughout the MVL trajectory until the saturation point, indicating a different spatial behavior of the members. This model achieves the largest dynamical spread among members, which is sustained until saturation. 5.2 MVL fingerprint: model spatiotemporal dynamics The dynamics of each single-model ensemble is better appreciated at the initial steps of the model simulation, while the initial perturbation is grown by the model. However, the DEMETER project aims at predicting climate features at longer timescales. At a global scale, the information contained in the initial conditions of the atmosphere is lost in the moment the MVL trajectory reaches saturation and perturbations become climatic fluctuations (the onset of saturation is approximately 15– 20 days). After that time, the MVL trajectory stays at a fixed point (for an stationary model) indicating the characteristic structure of climatic fluctuations. In the case of atmospheric models, the model is not stationary and there are changing forcings (the season, the slowly varying ocean coupled component, etc.) which determine a region where the trajectory ends up (see Fig. 7). The size and position of this region depends again on the model internal dynamics, but also on the forcings the model is subject to at these longer time scales. The analysis of the fluctuations in that region requires further research, since the stationary state arises from an equilibrium between the multiplicative growth induced by chaos and the saturation effect of
241
Fig. 7 Zoom over the saturation region (as outlined in Fig. 3) of the MVL diagram of the mean ensemble evolution for the November initializations. For clarity, the leading month (the entry path of each curve) is not plotted
non-linearities. When any of these effects varies slowly in time the state in the diagram moves accordingly. In any case, given that the forcings of each model are similar, the final region where the model ends is also an indication of the dynamics since it characterizes the amplitude and spatial structure of the climatological fluctuations of the model. Thus, models with similar dynamics will develop similar climatological fluctuations. This fact is empirically shown in Fig. 7, displaying a zoom over the region where the perturbations have saturated. There are four different groups of models: (1) the cnrm and crfc models, (2) the scnr model, (3) the smpi model and (4) the lody, scwf and ukmo. The two models in group (1) are the only models sharing the same atmospheric and similar oceanic components (see Table 1); two of the models of the main group (4) share the same atmospheric module, and the models in (2) and (3) have similar atmospheric component (different from the rest of models), but differ in the ocean. Thus, although these results are just empirical and need further interpretation and validation, they offer a promising tool for model intercomparison considering the different dynamical behaviors. The weighting of models to extract the most reliable information out of a multi-model ensemble is a challenging topic. The current approaches weight models according to the skill in reproducing certain variables (Krishnamurti et al. 1999; Yun et al. 2003; Pavan and Doblas-Reyes 2000; Doblas-Reyes et al. 2005; Stephenson et al. 2005) or, most commonly, an equal-weight approach is followed. There have been attempts to optimally combine models avoiding the co-linearity of the predictions by using Principal Component Analysis
123
242
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS
(Doblas-Reyes et al. 2005) or Maximum covariance analysis (Stephenson et al. 2005). Co-linearity is assessed in terms of model output variables. However, these methods focus on the skill of the different models to reproduce a particular observed variable, or a set of variables, and do not explicitly take into account that current multi-model ensembles are not composed of completely independent models. For instance, if the skill of all models forming DEMETER were similar, then double weight would be assigned to the ARPEGE, ECHAM and IFS atmospheric models with respect to the HadAM model, thus biassing the sampling process used to form an ensemble. The main idea behind multi-model ensemble prediction is exploring uncertainties in the dynamics, but using the same model, even in different versions, does not seem to contribute to the diversity in the dynamics and to fully explore the existing uncertainty (for instance, if only two models were to be included in the ensemble, maximum uncertainty would be attained combining models from different groups 1–4 described above). As we have shown in this paper, the MVL diagram is a promising tool to assess if a given GCM (using a given version of the atmospheric and oceanic components) adds dynamical diversity to the ensemble or just increases the ensemble balance towards a single-model ensemble (regardless of how the output variables compare to each other or with reality). Further research is needed to develop appropriate statistical methods to incorporate this information as a pre-screening process to be applied before any weighting scheme, in order to preserve diversity and uncertainty.
6 Summary and conclusions We applied systematically the MVL representation of ensemble perturbation growth to the model simulations carried out within the EU-funded DEMETER project. This representation unveils dynamical and ensemble generation features of the simulations whilst it is mostly insensitive to different model realizations (for different years or seasons). As a result, the model ensemble generation procedure can be tested in terms of its consistency with the model attractor and the assimilated or non-assimilated character of the initial perturbed conditions. Also, the dynamical similarity among different modeling systems can be assessed. A great part of the present work focuses on the initial stages (2–3 weeks) of the DEMETER model simulations. For the purpose of seasonal forecasting, these initial stages are completely disregarded, since seasonal forecasting focuses on larger time-scales and the uncertainties in the initial stages are much better sampled in short and medium range forecasting systems. We used the results from
123
DEMETER since it is an unprecedented and systematic exploration of the amplification of initial perturbations by a wide range of state-of-the-art GCMs. Thus, our results referring to the initial stages cannot be interpreted as a criticism to the project design. Even though there appears to be a strong impact of the imposed atmospheric conditions on the first few forecast days, the differences are not relevant after this 2–3 weeks period, which anyway was not regarded in the validation process. Our results concerning the model climatological state are, however, fully relevant for DEMETER. The similarities in the building blocks of the different models compromise the dynamical diversity of the ensemble, which was one of the main goals pursued by the project. Combining the models with an equal weight would lead to favoured building blocks. The MVL diagram presented in this paper and, in particular, the final saturation region provides a qualitative tool to assess the dynamical diversity among different models. The task to derive a quantitative measure for an optimal combination of models from the point of view of the dynamics is out of the scope of the present paper, but it is a promising field to explore in the future. Acknowledgments The authors wish to thank F.J. Doblas-Reyes for his helpful comments. The authors made extensive use of open source software for processing and plotting purposes. Namely, we used the Climate Data Operators (CDO) developed at the Max Planck Institute for Meteorology in Hamburg and the Generic Mapping Tools (GMT) developed at the University of Hawaii (Wessel and Smith 1991). This work was partially supported by the Spanish Ministry of Education and Science through grants CGL-2007-64387/CLI and CGL200506966-C07-02/CLI. J. F. is supported by the Spanish Ministry of Education and Science through the Juan de la Cierva program.
References Doblas-Reyes FJ, Hagedorn R, Palmer T (2005) The rationale behind the succes of multi-model ensembles in seasonal forecasting—II. Calibration and combination. Tellus 57A:234–252 Giorgi F, Mearns LO (2003) Probability of regional climate change based on the reliability ensemble averaging (REA) method. Geophys Res Lett 30(12):1629 Gneiting T, Raftery AE (2005) Weather forecasting with ensemble methods. Science 310:248–249 Gutie´rrez JM, Primo C, Rodrı´guez MA, Ferna´ndez J (2008) Spatiotemporal characterization of ensemble prediction systems - the mean-variance of the logarithms (MVL) diagram. Nonlin Processes Geophys 15(1):109–114. http://www.nonlinprocesses-geophys.net/15/109/2008 Hewitt CD, Griggs DJ (2004) Ensembles-based predictions of climate changes and their impacts. Eos Trans AGU 85(52):566. doi: 10.1029/2004EO520005 Krishnamurti TN, Kishtawal CM, LaRow TE, Bachiochi DR, Zhang Z, Williford CE, Gadgi S, Surendran S (1999) Improved weather and seasonal climate forecats from multimodel superensemble. Science 285:1548–1550 Lo´pez JM, Primo C, Rodrı´guez MA, Szendro I (2004) Scaling properties of growing noninfinitesimal perturbations in spacetime chaos. Phys Rev E 70:056224
J. Ferna´ndez et al.: MVL spatiotemporal analysis for model intercomparison in EPS Molteni F, Buizza R, Palmer T, Petroliagis T (1996) The ECMWF ensemble prediction system: methodology and validation. Quart J Roy Meteor Soc 122:73–119 Palmer TN (2002) The economic value of ensemble forecasts as a tool for risk assessment: from days to decades. Quart J Roy Meteor Soc 128:747–774 Palmer TN, Shukla J (2000) Editorial to DSP/PROVOST special issue. Quart J Roy Meteor Soc 126:1989–1990 Palmer TN, Alessandri A, Andersen U, Cantelaube P, Daveyand M, De´le´cluse P, De´que´ M, Dı´ez E, Doblas-Reyes FJ, Feddersen H, Graham R, Gualdi S, Gue´re´my JF, Hagedorn R, Hoshen M, Keenlyside N, Latif M, Lazar A, Maisonnave E, Marletto V, Orfila APMB, Rofel P, Terres JM, Thomson MC (2004) Development of a European multimodel ensemble system for seasonal-to-interannual prediction DEMETER. Bull Am Meteor Soc 85:853–872 Pavan V, Doblas-Reyes FJ (2000) Multi-model seasonal hindcasts over the Euro-Atlantic: skill scores and dynamic features. Climate Dyn 16:611–625 Primo C, Szendro I, Rodrı´guez MA, Gutie´rrez JM (2007) Error growth analysis in systems with spatial chaos: coupled map lattices and global weather models. Phys Rev Let 98:108501 Primo C, Rodrı´guez MA, Gutie´rrez JM (2008) Logarithmic bred vectors. A new ensemble method with adjustable spread and calibration time. J Geophys Res 113:D05116. doi:10.1029/ 2007JD008998
243
Raftery AE, Gneiting T, Balabdaoui F, Polakowski M (2005) Using bayesian model averaging to calibrate forecast ensembles. Mon Wea Rev 133:1155–1174 Stephenson DB, a SCoelho C, Doblas-Reyes FJ, Balmaseda M (2005) Forecast asssimilation: a unified framework for the combination of multi-model weather and climate predictions. Tellus 57A:253–264 Toth Z, Kalnay E (1993) Ensemble forecasting at NMC: the generation of perturbations. Bull Am Meteor Soc 74:2317–2330 Uppala SM, Ka˚llberg PW, Simmons AJ, Andrae U, Bechtold VDC, Fiorino M, Gibson JK, Haseler J, Herna´ndez A, Kelly GA, Li X, Onogi K, Saarinen S, Sokka N, Allan RP, Andersson E, Arpe K, Balmaseda MA, Beljaars ACM, Berg LVD, Bidlot J, Bormann N, Caires S, Chevallier F, Dethof A, Dragosavac M, Fisher M, Fuentes M, Hagemann S, Ho´lm E, Hoskins BJ, Isaksen L, Janssen PAEM, Jenne R, Mcnally AP, Mahfouf JF, Morcrette JJ, Rayner NA, Saunders RW, Simon P, Sterl A, Trenberth KE, Untch A, Vasiljevic D, Viterbo P, Woollen J (2005) The ERA-40 re-analysis. Quart J Roy Meteor Soc 131:2961–3012 Wessel P, Smith WHF (1991) Free software helps map and display data. EOS trans AGU 72:441 Yun WT, Stefanova L, Krishnamurti TN (2003) Improvement of the multimodel superensemble technique for seasonal forecasts. J Clim 16:3834–3840
123
