Multi-objective environmental model evaluation by means of multidimensional kernel density estimators: Efficient and multi-core implementations

1 This is the accepted manuscript of the article that appeared in final form in Environmental Modelling & Software 63 : 123-136 (2015), which has been published in final form at https://doi.org/10.1016/j.envsoft.2014.09.019. © 2014 Elsevier under CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/) Multi-objective environmental model evaluation by means of multidimensional kernel density estimators: efficient and multi-core implementations Unai Lopez-Novoa1, Jon Sáenz2,3, Alexander Mendiburu1, Jose Miguel-Alonso1, Iñigo Errasti4, Ganix Esnaola2, Agustín Ezcurra2, Gabriel Ibarra-Berastegi4 Corresponding author: Unai Lopez-Novoa {E-mail: [email protected], Phone: +34 943 018 012} 1Intelligent Systems Group, Dept. of Computer Architecture and Technology, University of the Basque Country (UPV/EHU), P. Manuel Lardizabal 1, 20018, Donostia-San Sebastian, Spain 2Dept. Applied Physics II, University of the Basque Country (UPV/EHU), Sarriena Auzoa z/g, 48940Leioa, Spain. 3Plentzia Marine Station (PIE-UPV/EHU), Areatza Pasealekua, 48620, Plentzia, Spain 4Dept. Nuclear Engineering and Fluid Mechanics, University of the Basque Country (UPV/EHU), Alda. Urkijo, s/n, 48013, Bilbao, Spain Abstract We propose an extension to multiple dimensions of the univariate index of agreement between PDFs used in climate studies. We also provide a set of high-performance programs targeted both to single and multi-core processors. They compute multivariate PDFs by means of kernels, the optimal bandwidth using smoothed bootstrap and the index of agreement between multidimensional PDFs. Their use is illustrated with two case-studies. The first one assesses the ability of seven global climate models to reproduce the seasonal cycle of zonally averaged temperature. The second case study analyzes the ability of an oceanic reanalysis to reproduce global sea surface temperature and sea surface height. Results show that the proposed methodology is robust to variations in the optimal bandwidth used. The technique is able to process multivariate datasets corresponding to different physical dimensions. The methodology is very sensitive to the existence of a bias in the model with respect to observations. Keywords: Multivariate kernel density estimation, multidimensional kernel density estimation, multi-core implementation, environmental model evaluation 2 Software availability: Name or software or dataset: density-parallel Developers: Unai Lopez-Novoa (1) Jon Sáenz (2) Alexander Mendiburu (1) Jose Miguel-Alonso (1) (1) Intelligent Systems Group, Dept. of Computer Architecture and Technology, University of the Basque Country (UPV/EHU), P. Manuel Lardizabal 1, 20018, Donostia-San Sebastian, Spain. [email protected], Phone: +34 943 018 012 , [email protected], Phone: +34 943015020 and [email protected] Phone: +34 943 018019 (2) EOLO Group, Department of Applied Physics II, University of the Basque Country, (UPV/EHU), Barrio Sarriena s/n, 48940-Leioa, Spain. [email protected], Phone: +34 946012445. Year first available: 2014 Hardware required: any single-core or multi-core processor. Software required: •C Compiler (e.g. GCC). OpenMP compiling capabilities required for multi-threaded implementations •MESCHACH math library. Available at http://homepage.math.uiowa.edu/~dstewart/meschach/ •netCDF data handling library. Available at http://www.unidata.ucar.edu/software/netcdf/ Program language: ANSI-C Program size, including example data: 102 KB in a single .tar.gz file. License: Open software, available under a “New BSD 3-clause license”. Availability: http://www.sc.ehu.es/ccwbayes/isg/index.php?option=com_remository&Itemid=13 Or go to “Downloads” section of http://www.sc.ehu.es/isg 3 1 Introduction Climate models are the best tools that scientists currently have in order to assess the impact of increasing concentrations of greenhouse gases and other anthropogenic influences on the observed climate of the Earth. These tools allow scientists to understand climatic changes from a dynamical point of view and to give quantitative answers to questions about future climate. They make feasible the assessment of the characteristics (deterministic versus stochastic) of some climatic variations at different temporal or spatial scales. Finally, they are fundamental in the attribution phase of the study of the climate change problem, since they allow to confidently discard competing hypothesis such as whether the climate change is rooted in natural or anthropogenic causes (Bengtsson, 2013; Knutti, 2008). Climate models must be evaluated against different observations (Otto et al., 2013) or paleoclimate data (Braconnot et al., 2012; Hind et al., 2012; Moberg, 2013; Sundberg et al., 2012) in order to get a quantitative indication on the confidence that we can put into their outputs depending on the efficacy of the models to reliably represent the climatic processes and feedbacks. Contrary to operational weather forecast models, there is no way to properly evaluate models against future climate, since future climate does not exist yet (Randall et al., 2007; Stocker et al., 2013). Additionally, since parameterizations of sub-grid scale processes are not fully independent from current climate, it is clear that evaluation against current climate is the only feasible, albeit not perfect, solution in terms of evaluating the projections for future climate (Errasti et al., 2011, 2013; Radić and Clarke, 2011; Reichler and Kim, 2008). The rationale behind this hypothesis is that models that are able to better simulate current climate are the ones that we expect will also be the best ones in terms of the simulation of future climate. It is well known that this is not necessarily true due to the different behavior of models in terms of their internal feedback mechanisms (Andrews et al., 2012; Dessler, 2013). These feedbacks lead to differing values of the climate sensitivity of models and, hence, future warming is dependent on these different sensitivities, leading to the question whether all the models are equally valid (Knutti, 2010). Particularly for downscaling applications and regional impact analysis, an evaluation of the adequacy of models is a common step (Brands et al., 2011; Radić and Clarcke, 2011; Walsh et al., 2008). Considering that numerical downscaling is computationally expensive, it cannot be 4 performed over all the models available from a big experiment such as CMIP3 or CMIP5, and it is usually only performed on a subset of the models (Hewitt and Griggs, 2004). Then, the best models are only used for downscaling in regional climate assessment, since they have been proven to be the most suitable over a particular region. This strategy of selecting the best subset of all the available models has been contested by some studies (Reifen and Toumi, 2009) and defended by others (Macadam et al., 2010), since this result depends on the removal of the seasonal cycle and the use of anomalies instead of raw output from the models. There are several inter-comparison experiments that have been set-up in order to drive the models with common boundary conditions so that results between model runs can be compared. As an example of this kind of standard experimental setups, that in several cases have their origins in the nineties, we can cite the Atmospheric Model Intercomparison Project (Gates et al., 1998), the Project for Intercomparison of Land Surface Parameterization Schemes (PILPS) (Henderson-Sellers et al., 1995) or the Palaeoclimate Modelling Intercomparison Project (PMIP) (Kageyama et al., 1999), the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP), described by Lamarque et al. (2013), or the one that is most relevant for our study, since we use data from this experiment, the Coupled Model Intercomparison Project (CMIP) (Meehl et al., 2007; Taylor et al., 2012). These projects have covered several phases through the years and for the case of the CMIP data, CMIP5 can already be used. In general, models grouped under a similar experimental set-up such as CMIP3 or CMIP5 are considered as an ensemble of opportunity (Annan and Hargreaves, 2010). There are some limitations because even for coordinated experiments such as CMIP3, there is some freedom in the way the external boundary conditions are applied (they are not 100% equal for all the models), see Table 1 in Wang et al. (2007). The number of realizations from every model in the ensemble of opportunity is not the same, neither and, therefore, the influence of each model in the behaviour of the ensemble is not the same. Additionally, models are less independent than they should be, since critical algorithms or components are shared by several models (Fernández et al., 2009; Knutti et al., 2013; Masson and Knutti, 2011; Pennell and Reichler, 2011). In terms of evaluation of climate models, it is well known that climate simulations are run, most of the times, past the limit of deterministic predictability associated with predictability of the first kind, according to Lorenz's classification. Climate models simulate climate change under varying boundary conditions in terms of the Probability Density Functions (PDF) of climatic variables. The varying boundary conditions consist of external forcings such as the variability in solar irradiance, 5 orbital parameters or anthropogenic greenhouse gas emissions, amongst other potential driving factors (Bengtsson, 2013). Some aspects of climate simulations are deterministic, such as the seasonal cycle at extratropical latitudes (Errasti et al., 2013). On the other side, some properties of the atmospheric circulation such as the blocking at extratropical latitudes can not be precisely forecast with lead times corresponding to several days without the use of ensemble forecast systems (Marshall et al., 2014) because of the sensitivity to initial conditions (Frederiksen et al., 2004) and the model formulation (Pelly and Hoskins, 2003) too. Therefore, climate model evaluations that are run past the limit of deterministic predictability should not a priori expect a close consistency between weather-linked variations of global or regional temperature or precipitation between models and observations, except at the longer time-scales responding to external forcings (Gleckler et al., 2008; Santer et al., 2011). These differences reflect the well-known difference of the sensitivity of the results to errors in the initial conditions (predictability of the first kind) or to errors in the evolving boundary conditions (predictability of the second kind) (Chu, 1999; Lorenz, 2006). There is not a universally accepted strategy for climate model evaluation, since it is well known that climate model evaluation and corresponding skill scores fundamentally depend on the target area, variable or intended application of the model evaluation study (Knutti, 2008). Some studies make a focus on the deterministic parts of the model simulations (basically, the seasonal cycle) (Boer and Lambert, 2001; Taylor, 2001). Other studies (oriented to the study of droughts or floods) tend to focus on extreme percentiles, since they are much more meaningful indicators of climate change (DeAngelis et al., 2013). This study by DeAngelis et al (2013) is currently important for us, because it shows that models sometimes produce accurate values for the average of some climatic variable due to error compensation effects. They can, for instance, underestimate the frequency of high precipitation events and overestimate the frequency of low precipitation events. This points to the need to evaluate additional characteristics of climate model simulations beyond the mean value and standard deviation. Consequently, a few years ago, an index computed from the whole PDF of climatic variables was developed (Maxino et al., 2008; Perkins, 2007). It compares two PDFs and computes the minimum value of both PDFs at every abscissa. The area below this minimum represents the area below both PDFs. As such, for a perfect model, its value would be one, if both PDFs matched perfectly. This PDF-index or index of distributional agreement is the one that we will generalize to the multidimensional case in this contribution. The PDF-index analyses the 6 correspondence of the whole PDF both from a model and observations (Maxino et al., 2008; Perkins et al., 2007). The one-dimensional PDF-index has very often been used in the literature through the last years (Brands et al., 2012; Errasti et al., 2011, 2013; Fu et al., 2013; Maxino et al., 2008; Perkins et al., 2007; Schwalm et al., 2013; Ylhaïsi and Räisänen, 2013 to name a few). In this contribution we propose its extension to multiple dimensions, thus allowing to compare several features of climate or environmental models at a single step. The use of the PDF-index shows some advantages with respect to other approaches, in the sense that it samples the full PDF of the climatic variables. Therefore, the PDF-index is a very good index for the overall evaluation of the agreement between climate models and observed climate. However, there are other shortcomings, such as the fact that the analysis in terms of PDFs does not consider the time sequence of events, and the number of frost days or the number of continuous days without precipitation are important in terms of impacts (Brands et al., 2012). The PDF-index gives less weight to the tails of the distribution and, it is thus not adequate as the single index for the analysis of extremes (Brands et al., 2012). In the case of the papers mentioned previously, the PDF-index is computed by means of unidimensional PDFs. However, in several cases, studies using the PDF-index or other scores (Dessai et al., 2005; Reichler and Kim, 2008) evaluate the skill of climate models according to several variables that may be of interest for the impact community. The most obvious instances might be precipitation and temperature, but if the scientists are interested in downscaling strategies, other variables such as geopotential height or sea level pressure appear very often (Brands et al., 2011; Errasti et al., 2011, 2013; Fu et al., 2013; Maxino et al., 2008; Radic and Clarke, 2011). In these previous references, the skill of the models is computed on a per-variable basis by means of univariate diagnostics. Their final skill score is computed by aggregating individual per-variable evaluations either by simple averaging or ranking of skill scores. However, there is currently a lack of universally accepted way of performing this combination of scores for different variables and the methodology that we propose in this contribution is aimed to fill this void, since a single index of distributional agreement is returned from the multidimensional PDF. The main objective of this contribution is, therefore, to develop a methodology that can be applied to get a multidimensional score that allows to evaluate in a single step different variables from climate simulations against observations. In order to explain the advantages derived from using a multidimensional approach, we show a simple example derived from a synthetic dataset. We have created three synthetic datasets, G1, G2 and G3, derived from two-dimensional gaussian 7 distributions. For each case, the gaussians are centered µ i = 0 but the corresponding covariance matrices used to create them are given by S 1 =S 2 = ( 1 0.750.75 1 ) for G1 and G2, whilst for G3, the third gaussian, the covariance matrix is given by S 3 = ( 1 −0.75−0.75 1 ) . It can be seen that, despite the difference in the structure of the distributions of points (Figure 1, left), the univariate PDFs and corresponding indices show a good agreement (Figure 1, middle and right and Table 1), even though the distributions are different. This is quite an artificial example, but it illustrates the point that some parts of the PDFs close to the diagonals can be projected onto similar areas over the axis when using unidimensional indexes of agreement, masking the differences between the PDFs of the model and the observations. Therefore, it is interesting to analyse the full structure of the multidimensional PDF, since it yields a realistic difference between the score corresponding to G1 versus G2 (good agreement) and G1 versus G3 (bad agreement). Figure 1. Points created from G1 (red) and G3 (green) distributions (left), univariate probability distributions corresponding to the X variables (middle) from G1 (red) and G3 (green) and univariate probability distributions corresponding to the Y variable from G1 (red) and G3 (green). Table 1. Indices of distributional agreements for points derived from known gaussians using univariate scores for X and Y variables or two-dimensional scores. Univariate score 2D score G1-X G1-Y G1 G2 X 0.998 0.999 Y 0.999 G3 X 0.998 0.463 Y 0.997 8 To fill in the gap illustrated by the example, we generalize the PDF-index by Perkins et al. (2007) to n-dimensional phase spaces with the final aim of allowing an easy multi-criteria evaluation of models. That way, a single PDF-based index can group the performance of the models according to the multidimensional phase-space spanned by all the variables chosen for the evaluation of the model. In order to make it easier for other researchers to use this methodology, we present an implementation of this multidimensional extension by means of a set of tools that can properly address the computational problems that appear when making kernel-based estimations of PDFs with massive datasets. In the case studies that we show in this paper, we compute several realizations of the PDF-index using one to four dimensional phase spaces with up to 13 000 points for every particular model/realization. This is very intensive computationally, and for this reason we feel that an efficient implementation of the estimation of multidimensional PDFs could be of great help for researchers in this area. Thus, we have developed a general-purpose tool to compute kernel-based multidimensional PDF estimations that runs on state-of-the-art multi-core processors. Our proposal has two main characteristics: (1) a fine-tuned algorithm to calculate the PDF that minimizes the number of computations and (2) a parallel implementation of this algorithm that allows it to efficiently run in multi-core processors. The method is applied to two different case-studies. The first case study corresponds to a realistic application of climate model evaluation (Errasti et al., 2013). Their results are re-analyzed using this new methodology. Additionally, a sensitivity of the results to the selection of the bandwidth parameter is carried out. Finally, the robustness of the results provided by the method to the existence of biases in the models is also studied. The second case study corresponds to the analysis of the performance of a coupled atmosphere-ocean reanalysis in reproducing the global scale Sea Surface Temperature and Sea Surface Height and it corresponds to a higher-dimensional problem. This second case study will be used to stress two of the merits attributed to the proposed methodology. On the one hand, this example in the context of the physical oceanography will demonstrate the wide range of the applicability of the methodology. On the other, the combination of variables with different physical dimensions will illustrate the ability of the method to process multivariate data and its ability to be applied to a large family of environmental models. The remaining of this paper is structured as follows. Section 2 presents the materials and methods used in the paper. Results are shown in section 3. The discussion is presented in section 4, and the paper finishes with conclusions in section 5. 9 2 Material and methods 2.1 Data representing the daily seasonal cycle of zonally averaged temperature from global climate models and reanalysis. For the first case study shown in this paper, we select a reduced dimensionality representation of the daily seasonal cycle of temperature that we already analyzed in Errasti et al. (2013). We analyze temperature of the air at the surface (TAS) daily data from seven models (20C3M simulations) of the CMIP3 experiment (Meehl et al., 2007) that were used by the Fourth Assessment Report of the IPCC (Randall et al., 2007). The models used are the BCCR-BCM2.0, GFDL-CM2.0, GFDLCM2.1, MIROC3.2-HR, MIROC3.2-MR, MPI-ECHAM5 and MRI-CGCM2.3, and the basis for the selection of this subset of models and their characteristics can be found in Errasti et al. (2013). The same procedure is used for TAS data from ERA40 (Uppala et al., 2005) and NCEP/NCAR Reanalysis 1 (Kalnay et al., 2996), referred to as NCEP onwards. The TAS data from models and reanalyses were re-gridded to the same 2.5º x 2.5º grid by means of bilinear interpolation, since this was the coarsest grid used by any of the reanalysis used (the one used by NCEP). To reduce the dimensionality of such a gridded dataset, the TAS data were zonally averaged and projected onto Legendre polynomials that constitute an adequate basis over the sphere. Those have very often been used as a basis in one-dimensional energy balance models (North et al., 1981) and are the basis used for the meridional components of spherical harmonics used in spectral decompositions of the equations of motion (Washington and Parkinson, 2005). The Legendre polynomials are orthonormal in the set of functions over the continuous interval [-1, 1], but their corresponding discrete-grid counterparts are not orthogonal. For this reason, we applied the Gram-Schmidt orthogonalization procedure to obtain the leading discrete orthogonal  P 0 ( µ ) ,  P 1 ( µ ) and  P 2 ( µ ) Legendre polynomials. In the previous equations, µ=sin ( θ ) refers to the sine of latitude. The zonally averaged TAS profiles have been projected onto the orthogonal discrete  P 0 ( µ ) ,  P 1 ( µ ) and  P 2 ( µ ) Legendre polynomials and this has provided us with the corresponding time-varying coefficients c 0 ( t ) , c 1 ( t ) and c 2 ( t ) . Due to the meridional shape of the orthogonal discrete polynomials shown in Figure 1 in Errasti et al., (2013), the physical meaning of the temporal expansion coefficients can be easily understood. The coefficient c 0 ( t ) describes the seasonal evolution of global-mean temperature linked to the different distances from the earth to the Sun corresponding to the apogee or the perigee positions, c 1 ( t ) describes the seasonal evolution of summer-winter from one Hemisphere to the other and c 2 ( t ) describes the TAS differences between 16 In the first case study used in this paper the phase-space is tridimensional, but the implemented program allows the user to go up to any dimension (starting from one), as shown by the second case study used in this paper. The one-dimensional case is also covered by our program even though there are other implementations that can cover the one dimensional case. The novelty of our contribution lays on the generalization of the one-dimensional score to higher dimensions. In addition, standard OpenMP programming directives (Dagum and Menon, 1998) have been also included with the aim of exploiting the multi-core capability of present computers. The set of observed points will be equally distributed amongst the processors for their computation in parallel. This way, the workload is split among the available processors, reducing the execution time (almost) linearly to the number of cores used. 2.3.2 Selection of the optimal bandwidth The second step that we describe (although it should be the first step in the application of the programs) is to find the optimal bandwidth value for the PDF computed in the first step. It is well known that the computation of the optimal bandwidth to be used in PDF estimations using kernels is a critical step in obtaining reliable PDFs. There are two major strategies for the determination of the optimal bandwidth (Scott, 1992; Slverman, 1986): cross-validation (Duong and Hazelton, 2005) and smoothed bootstrap (Faraway and Jhun, 1990). The cross-validation approach leads to the convolution of the kernel with itself, a very tough mathematical problem for the Epanechnikov kernel with an open number of dimensions. It is usually solved by means of gaussian multiplicative kernels (Duong and Hazelton, 2005), but this wouldn't allow us to use the OPB approach explained in section 2.3.1 above. In our case, we have selected the use of smoothed bootstrap estimates of the optimal bandwidth, since it simplifies the generalization of the solution to an open number of dimensions in the multidimensional case for the non-multiplicative Epanechnikov kernel we are using. In order to produce the new estimations in the multidimensional case we take advantage of the fact that the kernel is spherically symmetric in the space corresponding to the spherically symmetric principal components. Thus, the same strategy used by univariate kernel estimations is used for every direction in the space spanned by the spherically symmetric principal components (Silverman, 1986). Surrogate samples are created in this space, and this procedure guarantees that the structure of the covariance matrix is properly preserved. 17 Following Faraway and Jhun (1990), we use a smoothed bootstrap procedure to estimate the squared error between two estimates of the PDF. The smoothed estimate starts from a reference evaluation of the PDF  f ( x,h 0 ) computed using a reference bandwidth h 0 . Then, several estimations  f n ( x,h ) of the PDF are performed at varying values of the bandwidth parameter h and a number of n= 1 … N realizations for every h. The bootstrap program checks the error between the “reference” PDF used in the smoothed bootstrap procedure and the actual bootstrap samples by evaluating the squared error ε n ( h ) = ∫ (  f ( x,h 0 ) −  f n ( x,h ) ) 2 dx . Then, the bootstrap-derived distribution of the squared errors is used to infer minimum, maximum, median, P 0 . 025 (2.5%) and P 0 . 975 (97.5%) percentiles of squared error for every value of h. This information is reported to the user at every h value. The h value producing the lowest values of the error estimates (we use the median of ε n ( h ) in the case studies in this paper) is the one selected as the optimum bandwidth. This procedure has been implemented in the mpdfestimator_bootstrap program. It takes as input all the observed points and optionally (1) a reference bandwidth value, (2) a range of bandwidth values to be evaluated, (3) the boundaries of the evaluation space, and (4) the number of repetitions for the random sampling. If (1) is missing, the default corresponding to a multidimensional gaussian distribution with the same sample size is applied. If (2) is missing a range (+/- 20% around (1), with a step such that the maximum bandwidth interval is divided in 10 subintervals) is defined. In case (3) is missing, mpdfestimator_bootstrap defines a range that ensures a space that surrounds all the observed points. Finally, if (4) is missing, 500 realizations are performed. The program generates as output squared errors for each of the provided bandwidth values. The pseudo-code is shown in Listing 2. Compute the PDF for the reference bandwidth h 0 for each h in the range [hmin,hmax]{ for iter=1 to max_repetitions{ generate a random sub-sample S compute PDF for S compute squared error } generate statistics of squared error } return statistics 18 Listing 2. Pseudo-code for the mpdfestimator_bootstrap program. 2.3.3 Computation of the PDF score The final step of the methodology is to compute the PDF score. Once the user has computed the PDFs (by means of the mpdfestimator program) corresponding both to the model and the observations using the optimal bandwidth value reported by mpdfestimator_bootstrap, program mpdf_score has to be executed to get the PDF score against the reference model. The program mpdf_score takes as input two n-dimensional PDFs stored as netCDF files, generated for the same domain by the first program mpdfestimator, and provides as output a PDF-index S by means of the following equation (adapted in this case for a three-dimensional example): S= ∑ min ( Z ijk o ,Z ijk m ) dx i dx j dx k , where Z ijk o and Z ijk m refer to the evaluation of the PDF from observations and the model, respectively. Please note that for higher dimensions, the extension is straightforward. The equation closely follows the one used by Perkins et al. (2007) or Maxino et al. (2008), but has been extended in this case for its use with a PDF defined in a multi-dimensional (n-dimensional) space. Additionally, when working in several dimensions, the volume of the n-dimensional interval where the PDF is being computed must be taken into account for normalization purposes, and so, the dx i , dx j and dx k terms account for the fact that the range of the different variables can be very different (the average standard deviations of the coefficients in our first case study are 1.9 K, 9.7 K and 2.1 K). The program warns the user in case the bias for any of the dimensions is greater than 5% of the standard deviation of that variate. 2.4 Representation of marginalized PDFs for the interpretation of results. Finally, even though it is not part of the methodology we propose, in order to be able to identify the differences in the index corresponding to individual models and for illustration purposes of the results corresponding to the first case study, we have computed marginalized  f 2D ( c i ,c j ,h ) = ∫  f ( x,h ) dc k ,i≠j≠k two-dimensional PDFs and projected them onto the i-j C0C1, C0-C2 and C1-C2 planes, after marginalizing k axes C2, C1 and C0, respectively. This will allow us to show that using a single multidimensional score is better than using a set of unidimensional scores. We only present marginalized PDFs for the first case study in the paper and, for the second case study we just collect the aggregated values of the score in a table. 19 3 Results 3.1 Application to climate model simulation of the daily cycle of surface temperature Figure 3 shows the evolution of the median of the squared errors and the 95% confidence interval computed from the bootstrap analysis corresponding to the ERA40 data when the reference PDF is computed with two conservative estimates of bandwidth (h 0 =0.8 and h 0 '=0.67) against the bandwidth that would correspond to the same sample size for a gaussian PDF, 0.637. It can be seen that the bootstrap estimate suggests a slightly lower value (0.55-0.60) of the bandwidth parameter than the one that would correspond to a gaussian multidimensional PDF. As will be identified in the marginalized PDFs later, this is to be expected, since the zonally averaged surface temperature is very non-normal and periodic, so that several fine scale features of the PDF must be resolved, and they can only be properly resolved if the bandwidth is not very high. Therefore, in the following steps an optimum bandwidth of h=0.6 will be used unless otherwise explicitly stated. In order to test the sensitivity of the classification to different values of the bandwidth used, Table 2 presents the results of the multidimensional PDF scores for different values of the bandwidth parameter (every model is centered and checked against ERA40). For the optimum bandwidth (h=0.6), the best model available is the alternative reanalysis that is used in this study (the NCEP). This is something that we expected from the beginning, since both reanalyses are based on observations. This result supports the use of the method, since the method yields better results for alternative observation-based reanalyses. MIROC3.2-MR and HADGEM1 are the models that follow. Some of the model runs differ only on the initial conditions and most of them are grouped together, with the exception of HADGEM1. The interpretation of this result is that the variability of the index to the use of different initial conditions is very low, as should be expected. MIROC3.2HR, GFDL, ECHAM5 and BCM2 follow the previous models. The ranking finishes (for the subset of models and diagnostic variable used in this study) by the five random runs corresponding to the MRI model. All the runs corresponding to MRI are grouped, with low values of the score that do not mix with values corresponding to the rest of the models. It seems, therefore, that the intraensemble variance is in general (without the exception of MIROC3.2-MR and HADGEM1) smaller than the inter-model variance of the score. In general, the main characteristics of these results are robust even with changes in the bandwidth that span a -33% to a +33% interval from the optimum value found by means of bootstrap. The models that show the highest (lowest) performances with the optimum value of the bandwidth continue showing a similar performance for higher or lower values of the bandwidth. There are occasional excursions of a model to at most one alternative position up/down of the ranking, but, on the whole, models tend to maintain their relative ranks 20 even when the bandwidth is changed by a +/-33% relative change around the optimum value. Figure 3. Squared errors (median and 95% confidence interval as derived from the bootstrap estimates) between the randomly generated PDFs and the reference PDF (left, h 0 =0.8 and right, h 0 =0.67) used for the generation of the smoothed bootstrap. Figures 4, 5 and 6 show the plots of the marginalized PDFs for the case of the NCEP (contours) versus ERA40 (shaded), a model showing a high value of the score, MIROC-3.2-MR (contours) versus ERA40 (shaded) and a model with a lower score, such as MRI (contours) versus ERA40 (shaded). In order to show simple numbers in the plots and scales, values of the marginalized PDFs are multiplied by one thousand before plotting. Before computing the PDFs, the biases between every model and ERA40 have been removed by centering all the series. Table 2. Values of the multidimensional S score corresponding to different values of the bandwidth parameter and associated rankings that would correspond to the models, when compared with ERA40 reanalysis data. 21 Figure 4. Marginal PDFs of NCEP (contour) and ERA40 (shaded) projected onto the planes defined by the C0-C1 coefficients (left), C0-C2 coefficients (middle) and C1-C2 coefficients (right). Values of the PDF have been multiplied by 1000 in order to improve the representation of numbers. Figure 4, left, shows that on the C0-C1 plane, the PDF is clearly bimodal, as should be expected from a periodic deterministic signal such as the seasonal cycle of temperature. C0 represents the global average of surface temperature and C1 represents the difference in temperature between the Northern and Southern Hemispheres. The main clusters of the C0-C1 PDF appear aggregated around each Hemisphere's summer. NCEP values show a slightly warmer global temperature (C0) during Southern Hemisphere summer than the values shown by ERA40. Figure 4 (middle) shows that the amplitudes and phases of the mean global temperature (C0) and the equatorial bulge (C2) are similar in both reanalyses. On the C1-C2 plane, the marginalized PDF shows that the main difference between both reanalyses appears as a slightly higher difference of temperature between hemispheres (C1) in NCEP when the coefficient representing the equatorial bell (C2) is positive 22 (summer in the Northern Hemisphere). However, the PDFs generated by both reanalysis are extremely similar, as reflected in the high value of the S index between NCEP and ERA40 (0.82). This is something that we expected from the beginning, since they correspond to observational datasets. Figure 5. Marginal PDFs of MIROC3.2-MR (run 2, contours) and ERA40 (shaded) projected onto the planes defined by the C0-C1 coefficients (left), C0-C2 coefficients (middle) and C1-C2 coefficients (right). Values of the PDF have been multiplied by 1000 in order to improve the representation of numbers. Figure 5 corresponds to the marginal PDFs for MIROC3.2-MR model (second run), one of the best CMIP3 models according to the metric selected in this study. Over the C0-C1 plane (left), there is quite a good agreement between both PDFs, since both clearly represent the bimodal structure of the PDF. However, the differences between MIROC3.2-MR and ERA40 are higher than in the previous case, both in terms of the location of the Northern Hemisphere summer and also in transitions between seasons that appear in the areas between the maxima in the marginal PDF. In the case of the C0-C2 plane (middle), the highest disagreement appears at the precise location of the maxima of the marginal PDFs, particularly during Northern Hemisphere summer. A similar diagnostic can be derived from the marginal PDF over the C1-C2 plane. Despite both marginal PDFs are clearly bimodal, slight differences exist at the placing of the PDF maxima. The equatorial bell (C2) in MIROC3.2-MR is stronger than the one in ERA40 during negative phases (Northern Hemisphere winter) of inter-hemispheric temperature differences (C1). 23 Figure 6. Marginal PDFs of MRI (run 1, contours) and ERA40 (shaded) projected onto the planes defined by the C0-C1 coefficients (left), C0-C2 coefficients (middle) and C1-C2 coefficients (right). Values of the PDF have been multiplied by 1000 in order to improve the representation of numbers. Figure 6 corresponds to MRI (run 1) model. The evolution of the daily seasonal cycle of temperature in terms of C0 (global T) and C1 (inter-hemispheric temperature contrast, left) does not present a bimodal structure with the PDF maxima placed at the same points shown by the reanalysis. The transitions between summer and winter regimes happen through routes that do not correspond to the ones in the ERA40 Reanalysis. The structure of the marginal PDF for the C0-C2 plane is markedly different between MRI and ERA40, with the cold maximum in the PDF during summer in the Southern Hemisphere quite misplaced in the case of MRI. This is also apparent in the marginal PDF corresponding to the C1-C2 plane, where maxima of the PDFs do not appear neither on the same places nor even with the same phases. Finally, Figure 7 shows that the index is very sensitive to the existence of a bias between the models and reference observations. In this case, the PDFs are computed without previously removing the bias between the surrogate model (NCEP data) and the observations (ERA40) and the S score index that we get between ERA40 and NCEP reanalyses is extremely low (S=0.075). The marginal PDFs show that in general there is a very good agreement in the structure of the 3D PDFs, but the center of masses of both PDFs are not located at the same places. The biases for every coefficient are not very high, considering their variances. The bias of the C0 component is 0.7 K (0.2% relative error), the bias in C1 is 0.4 K (7.5% relative error) and the bias in C2 is -0.34 K (-1.33% relative error). However, even such low values of the bias lead to a score index that could be interpreted as poor performance of the surrogate model (NCEP reanalysis) versus ERA40 due to the complex structure of the 3D PDF. However, this is a false impression that can not be defended if the spatial patterns of 24 the marginal PDFs are analyzed in detail. This means that the index should not be applied to model results that are biased against the reference observations. The existence of biases in the models leads to greater observational uncertainty when the model and observational datasets are not centered. The code does not force the centering of the datasets and, therefore, the user must take care of this when the dimensionality reduction stage of the data analysis is done. In particular, it is interesting to stress that, internally, when computing the n-dimensional PDFs, all the datasets are centered (each one using its n-dimensional average) before computing the corresponding spherically symmetric principal components. When the output netCDF files holding the PDFs are saved, the original units in the phase space of each dataset (model or observations) are recovered and the average is added to the anomalies derived from the PDF in the principal component space stored in the memory of the computer. Therefore, the key point here is that if there exists a constant bias between the model and the reference observations (first and second netCDF files passed to program mpdf_score), it could lead to very low values of the score despite the model representing properly the variability (anomalies). This means that the evaluation of the models in terms of a constant bias and the n-dimensional PDFs should be carried out as different steps. Figure 7. Marginal PDFs of non-centered NCEP (contours) and ERA40 (shaded) projected onto the planes defined by the C0-C1 coefficients (left), C0-C2 coefficients (middle) and C1-C2 coefficients (right). Values of the PDF have been multiplied by 1000 in order to improve the representation of numbers. The bias between both reanalysis has been retained. From the point of view of performance, we have measured the execution time needed by each version of the program to complete the bootstrap procedure. On average, the serial OPB approach is 140 times faster than the serial GPB approach and, moreover, the parallel OPB program scales 25 linearly with the number of cores, being 4.3 times faster than its serial counterpart when using 4 cores. This means that an evaluation of a single model that takes approximately 22 days with the serial GPB program, can be executed in less than one hour using the most efficient and parallel implementation presented in this contribution. These experiments have been conducted in a desktop computer with an Intel i7 3820 processor (four cores, 3.6GHz, Hyperthreading enabled) with 8GB of RAM. Therefore, the use of this technique is not limited to the availability of specialized clusters or hardware that would limit its practical use. 3.2 Evaluation of Sea Surface Temperature and Sea Surface Height The first two PCs of the global coverage weekly time-scale SST (T1, T2) and SSH (H1, H2) variables belonging to the ARMOR-3D (Guinehut et al., 2004; Guinehut et al., 2012) and CFSR (Saha et al., 2010; Saha et al., 2014) datasets will be used in the following to evaluate the second with respect to the former. This means that the ARMOUR-3D product (blended satellite and in-situ observation product) is the reference to evaluate the CFSR product (coupled atmosphere-ocean modelling product). Considering the first two PCs of each variable in the evaluation (T1, T2, H1, H2), the global-scale main variability modes of each variable are being take into account at a glance. As the seasonal cycle was not explicitly removed from the anomalies used to deduce the PCs, the four considered variables are almost completely related to the global-scale seasonal cycle (H2 contains some longer time-scale variability). Thus comparing combinations of different variables from CFSR with those of ARMOR-3D the capacity of the modelling product to jointly characterize different main global-scale variability modes (their seasonal cycles) is evaluated. For example, if T1 and H1 are considered at the same time (case T1H1) the capacity of CFSR to simulate the main global-scale components of the seasonal cycle of the SST and SSH variables is being evaluated in a single and multivariate score. Table 3 shows the optimal h and the score obtained with the 6 analyzed cases going from the univariate T1 and H1 cases, the multi-dimensional univariate T1T2 and H1H2 (reserving the term multivariate to the cases with variables with different physical dimensions, i.e. Kelvins and meters) and the multivariate T1H1 and T1T2H1H2 cases. All variables have zero mean so no bias related issues will be observed in this case. Like in the previous case study on the TAS, the same three-step methodology was applied in this case: for a given row in Table 3, the optimal h using the bootstrap procedure is initially estimated. Next, the PDF using the optimal h is computed and, finally, the score (one dimensional, multidimensional or multivariate) is computed from the PDFs obtained from the CFSR and the ARMOR-3D variables. 32 The use of a multidimensional analysis produces a single index corresponding to every model even after analyzing several variables, and this result makes it easy to perform evaluation of the models under several target variables. In the contribution presented here, we have explored one case such as three Legendre coefficients that expand the daily cycle of zonally averaged temperature. However, the same approach could be applied to the joint analysis of temperature, outgoing longwave radiation or cloud cover (to name a few) such that the structure of the multidimensional PDFs (probably properly marginalized as in this contribution) could shed light over the behaviour of models according to known physical mechanisms. Even though a tool of the set described in this contribution allows to make an objective selection of the optimal bandwidth to be used in the generation of the PDFs by means of smoothed bootstrap, the case study in this paper shows that the ranking of the models is quite robust even under severe (+/- 30% of the optimal bandwidth) changes in the bandwidth used for the generation of the multidimensional PDFs. Thus, the results obtained through the use of the tools presented in this contribution are reliable. However, the index is extremely sensitive to the existence of a constant bias between models and it should not be used without previously centering the data, a finding in agreement with previous studies using univariate PDFs (Brands et al., 2011; Brands et al., 2012). The programs do not request that the datasets are centered, but a constant bias between the model and observations could lead to unphysical diagnostics in several dimensions. A potential solution is to perform the evaluation onto n-dimensional centered data, so that the bias is automatically removed. Alternatively, the analysis can be performed removing the bias from the model results with respect to observations. Therefore, the analysis of the bias must still be kept independent from the analysis of the shape of the PDF presented in this contribution. The current implementation of the mpdf_score program provides a warning if the bias at any of the dimensions is greater than 5% of the standard deviation of the correspoding variate. The second case study demonstrated the applicability of the proposed methodology to multivariate and multidimensional data using data from oceanographic SST and SSH variables too. In addition, and although it is not part of the proposed technique, this case study also demonstrated the potential of the use of a preprocessing step for the reduction of the dimensionality of the data, based on a PCA analysis of the original dataset in this case. This shows that the method can potentially be applied to a large family of environmental problems. 33 The overall evaluation of environmental models is a complex task and different performance scores detect different weak or strong points of the available global models. We hope that the addition of a new methodology and tools that allow its easy application by other researchers make it easier the identification in future experiments of areas of models that can be improved. Acknowledgements: Authors acknowledge constructive comments by three referees and the editor of this paper. These comments have lead to an improved version of the manuscript. We acknowledge the modeling groups, the Program for Climate Model Diagnosis and Inter-comparison (PCMDI) and the WCRP's Working Group on Coupled Modeling (WGCM) for their roles in making available the WCRP CMIP3 multi-model dataset. Support of this dataset is provided by the Office of Science, U.S. Department of Energy. ECMWF ERA-40 data used in this study have been provided by ECMWF. NCEP reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at http://www.esrl.noaa.gov/psd/ have been used. CFSR and CFSv2 data was provided by the Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, Boulder, Colorado. ARMOUR-3D data was obtainned from MyOcean (http://www.myocean.eu/). Authors thank financial funding by project CGL2013-45198-C2-1-R (MINECO, National R+D+i plan), the SAIOTEK program from the Basque Government (project S-P11UN137). Additional funding from different calls from the University of the Basque Country (UFI 11/55, PPM12/01 and GIU 11/01) has allowed this paper to be finished. This work has also been partially supported by the Saiotek and Research Groups 20132018 (IT-609-13) programs (Basque Government), TIN2010-14931 (Ministry of Science and Technology), COMBIOMED network in computational bio-medicine (Carlos III Health Institute). U. Lopez-Novoa holds a grant from the Basque Government. J. Miguel-Alonso and A. Mendiburu are members of the HiPEAC European Network of Excellence. Author contributions: JS, AM and JMA designed the research; ULN, JS, AM and JMA wrote the code distributed with this contribution; IE, AE, GIB and JS performed the computations that lead from data in the CMIP3 repository to the Legendre coefficients used in the first case study; GE performed the computations for the second example, ULN, JS, AM and JMA prepared the specific computations, graphics and results used in this contributions after the Legendre coefficients and ULN, JS, AM, IE, GE and JMA wrote the paper. References 34 Andrews, T., Gregory, J. M., Webb, M. J., and Taylor, K. E., 2012. Forcing, feedbacks and climate sensitivity in CMIP5 coupled atmosphere-ocean climate models. Geophysical Research Letters, 39, L09712, doi: 10.1029/2012GL051607. Ahamada, I. And Flachaire, E. 2010. Non-parametric Econometrics. Oxford University Press, 176 pages, Oxford. Annan, J. D., and Hargreaves, J. C. 2010. Reliability of the CMIP3 ensemble, Geophysical Research Letters, 37, L02703, doi:10.1029/2009GL041994. Bellman, R., 1961. Adaptive Control Processes: A Guided Tour . Princeton University Press, 255 pp. Bengtsson, L., 2013. What is the climate system able to do “on its own”? Tellus B65, 20189, doi:10.3402/tellusb.v65i0.20189. Bennett, N.D., Croke, B.F.W., Guariso, G., Guillaume, J.H.A., Hamilton, S.A., Jakeman, A.J., Marsili-Libelli, S., Newham, L.T.H., Norton, J.P., Perrin, C., Pierce, S.A., Robson, B., Seppelt, R., Voinov, A.A., Fath, B.D., Andreassian, V., 2013. Characterising performance of environmental models. Environmental Modelling & Software, 40, 1-20, doi:10.1016/j.envsoft.2012.09.011. Boer, G. J. and Lambert, S. J., 2001. Second-order space-time climate difference statistics. Climate Dynamics 17, 213-218, doi: 10.1007/PL00013735. Braconnot, P., Harrison, S. P. , Kageyama, M., Bartlein, P. J. , Masson-Delmotte, V., Abe-Ouchi, A., Otto-Bliesner, B., Zhao, Y., 2012. Evaluation of climate models using palaeoclimatic data , Nature Clim. Change , 2, 417-424, doi: 10.1038/nclimate1456 Brands, S., Herrera, S., San-Martín, D., Gutiérrez, J., 2011. Validation of the ENSEMBLES global climate models over southwestern Europe using probability density functions, from a downscaling perspective. Climate Research, 48, 145-161, doi: 10.3354/cr00995. Brands, S. Gutiérrez, J.M., Herrera, S., Cofiño, A. S., 2012. On the use of reanalysis data for 35 downscaling. Journal of Climate 25, 2517-2526, doi: 10.1175/JCLI-D-11-00251.1 Chu, Peter C., 1999. Two Kinds of Predictability in the Lorenz System. Journal of Atmospheric Sciences, 56, 1427–1432. doi: 10.1175/1520-0469(1999)056<1427:TKOPIT>2.0.CO;2 Dagum, L.; Menon, R., 1998 OpenMP: an industry standard API for shared-memory programming, Computational Science & Engineering, IEEE , 5, 46-55, doi: 10.1109/99.660313 DeAngelis, A. M., Broccoli, A. J., Decker, S. G., 2013. A Comparison of CMIP3 simulations of precipitation over North America with observations: Daily statistics and circulation features accompanying extreme events, Journal of Climate, 26, 3209-3230, doi: 10.1175/JCLI-D-12-00374.1 Dessai, S., Lu, X. and Hulme, M., 2005. Limited sensitivity analysis of regional climate change probabilities for the 21st century. Journal of Geophysical research, 110, D19108, doi: 10.1029/2005JD005919. Dessler, A. E., 2013. Observations of Climate Feedbacks over 2000–10 and Comparisons to Climate Models , Journal of Climate 26, 333-342, doi: 10.1175/JCLI-D-11-00640.1 Duong, T. and Hazelton, M. L., 2005. Cross-validation bandwidth matrices for multivariate kernel density estimation, Scandinavian Journal of Statistics 32, 485-506, doi:10.1111/j.14679469.2005.00445.x Errasti, I., Ezcurra, A., Sáenz, J., Ibarra-Berastegi, G., 2011. Evaluation of IPCC AR4 models over the Iberian Peninsula, Theoretical and Applied Climatology, 103, 61-79, doi: 10.1007/s00704-0100282-y Errasti, I., Ezcurra, A., Sáenz, J., Ibarra-Berastegi, G., Zorita, E., 2013. Comparison of the main characteristics of the daily zonally averaged surface air temperature as represented by reanalysis and seven CMIP3 models, Theoretical and Applied Climatology, 114, 417-436, doi: 10.1007/s00704-013-0842-z Faraway, J. J., Jhun, M., 1990. Bootstrap choice of bandwidth for density estimation. Journal of the American Statistical Association, 85, 1119-1122 36 Fasano, G., Franceschini, A. 1987. A multidimensional version of the Kolmogorov–Smirnov test, Monthly Notices of the Royal Astronomical Society, 225, 155-170, doi: 10.1093/mnras/225.1.155. Fernández, J., Primo, C., Cofiño, A. S., Gutiérrez, J. M., Rodríguez, M. A. 2009. MVL spatiotemporal analysis for model intercomparison in EPS: application to the DEMETER multimodel ensemble. Climate Dynamics, 33, 233-243, doi: 10.1007/s00382-008-0456-9 Frederiksen, J. S., Collier, M. A., Watkins, A. B. 2004. Ensemble prediction of blocking regime transitions. Tellus, 56A, 485-500, url: http://www.tellusa.net/index.php/tellusa/article/view/14460. Fu, G., Liu, Z., Charles, S. P., Xu, Z., Yao, Z., 2013. A score-based method for assessing the performance of GCMs: A case study of southeastern Australia. Journal of Geophysical research, 118, 4145-4167, doi: 10.1002/jgrd.50269 Gates, W. L. , Boyle, J. S., Covey, C., Dease, C. G., Doutriaux, C. M., Drach, R. S., Fiorino, M., Gleckler, P.J., Hnilo, J. J., Marlais, S. M., Phillips, T. J., Potter, G. L., Santer, B. D., Sperber, K. R., Taylor, K. E., Williams, D. N. 1999. An Overview of the Results of the Atmospheric Model Intercomparison Project (AMIP I). Bull. Amer. Meteor. Soc., 80, 29–55, doi: 10.1175/1520-0477(1999)080<0029:AOOTRO>2.0.CO;2 Gleckler, P. J., Taylor, K. E., Doutriaux, C., 2008. Performance metrics for climate models. Journal of Geophysical Research 113, D22105, doi: 10.1029/2007JD008972 Guinehut S., Le Traon, P.-Y., Larnicol, G., Philipps, S., 2004. Combining Argo and remote-sensing data to estimate the ocean three-dimensional temperature fields - A first approach based on simulated observations. Journal of Marine Systems, 46 (1-4), 85-98, doi: 10.1016/j.jmarsys.2003.11.022 Guinehut S., Dhomps, A.-L., Larnicol, G., Le Traon, P.-Y., 2012. High resolution 3D temperature and salinity fields derived from in situ and satellite observations. Ocean Science, 8(5):845–857, doi:10.5194/os-8-845-2012 Henderson-Sellers, A., Pitman, A. J., Love, P. K., Irannejad, P., Chen, T. 1995. The project for 37 Intercomparison of land surface parameterisaton schemes PILPS) Phases 2 and 3. Bull. Amer. Meteor. Soc., 76, 489-503, doi: 10.1175/1520-0477(1995)076<0489:TPFIOL>2.0.CO;2 Hewitt, C. D. and Griggs, D. J. 2004. Ensembles-based predictions of climate changes and their impacts, Eos Transactions of the AGU, 85, 566–566, doi:10.1029/2004EO520005. Hind, A., Moberg, A. Sundberg, R. 2012. Statistical framework for evaluation of climate model simulations by u se of climate proxy data from the last millennium – Part 2: A pseudo-proxy study addressing the amplitude of solar forcing, Climate of the Past, 8 1355-1365, doi: 10.5194/cp-8-1355-2012. Justel, A., Peña, D., Zamar, R. 1997. A multivariate Kolmogorov-Smirnov test of goodness of fit. Statistics and Probability Letters, 35, 251-259, doi: 10.1016/S0167-7152(97)00020-5 Kageyama, M., Valdes, P. J., Ramstein, G., Hewitt, C., Wyputta, U. 1999. Northern Hemisphere Storm Tracks in Present Day and Last Glacial Maximum Climate Simulations: A Comparison of the European PMIP Models. Journal of Climate, 12, 742–760, doi: 10.1175/1520-0442(1999)012<0742:NHSTIP>2.0.CO;2 Kalnay, E. , Kanamitsu, M. , Kistler, R. , Collins, W. , Deaven, D. , Gandin, L. , Iredell, M. , Saha, S. , White, G. , Woollen, J. , Zhu, Y. , Leetmaa, A. , Reynolds, R. , Chelliah, M. , Ebisuzaki, W. , Higgins, W. , Janowiak, J. , Mo, K. C. , Ropelewski, C. , Wang, J. , Jenne, R., Joseph, D., 1996. The NCEP/NCAR 40-year reanalysis project, Bulletin of the American Meteorological Society, 77, 437470. doi: 10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2 Knutti, R., 2008. Should we believe model predictions of future climate change? Philosophical Transactions of the Royal Society A, 366, 4647-4664, doi: 10.1098/rsta.2008.0169. Knutti, R., 2010. The end of model democracy? An editorial comment, Climatic Change 102, 395404, doi: 10.1007/s10584-010-9800-2 Knutti, R., Masson, D., Gettelman, A., 2013. Climate model genealogy: Generation CMIP5 and how we got there. Geophysical Research Letters, 40, 1194-1199, doi:10.1002/grl.50256. 38 Lamarque, J.-F., Shindell, D. T., Josse, B., Young, P. J., Cionni, I., Eyring, V., Bergmann, D., Cameron-Smith, P., Collins, W. J., Doherty, R., Dalsoren, S., Faluvegi, G., Folberth, G., Ghan, S. J., Horowitz, L. W., Lee, Y. H., MacKenzie, I. A., Nagashima, T., Naik, V., Plummer, D., Righi, M., Rumbold, S. T., Schulz, M., Skeie, R. B., Stevenson, D. S., Strode, S., Sudo, K., Szopa, S., Voulgarakis, A., and Zeng, G. 2013. The Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP): overview and description of models, simulations and climate diagnostics, Geosci. Model Dev., 6, 179-206, doi:10.5194/gmd-6-179-2013. Lopes, R. H. C., Hobson, P. R., Reid, I. D. 2008. Computationally efficient algorithms for the twodimensional Kolmogorov–Smirnov test. Journal of Physics: Conference Series, 119, 042019, doi:10.1088/1742-6596/119/4/0420 Lorenz, E. N. 2006. Predictability, a problem partly solved, Chapter 3 in Palmer, T. and Hagedorn, R. (eds.) Predictability of Weather and Climate, Cambridge University Press, Cambridge, 702 pp. Macadam, I., Pitman, A. J., Whetton, P. H., Abramowitz, G., 2010. Ranking climate models by performance using actual values and anomalies: Implications for climate change impact assessments. Geophysical Research Letters 37, L16704, doi: 10.1029/2010GL043877. Marshall, A.G., Hudson, D., Hendon, H.H., Pook, M. Alves,O., Wheeler, M. 2013. Simulation and prediction of blocking in the Australian region and its influence on intra-seasonal rainfall in POAMA-2. Climate Dynamics doi: 10.1007/s00382-013-1974-7. Masson, D., Knutti, R., 2011. Climate model genealogy, Geophysical Research Letters, 38, L08703, doi: 10.1029/2011GL046864. Maxino, C. C., McAvaney, B. J., Pitman, A. J., Perkins, S. E., 2008. Ranking the AR4 climate models over the Murray-Darling Basin using simulated maximum temperature, minimum temperature and precipitation. International Journal of Climatology, 28, 1097-1112, doi: 10.1002/joc.1612. Meehl, G. A., Covey, C., Taylor, K. E., Delworth, T., Stouffer, R. J., Latif, M., McAvaney, B., Mitchell, J. F. B., 2007. The WCRP CMIP3 multimodel dataset: A new era in climate change research. Bulletin of the American Meteorological Society, 88, 1383-1394. doi: 10.1175/BAMS- 39 88-9-1383 Menke, W. and Menke, J., 2012, Environmental Data Analysis with Matlab, 263 pages, Elsevier, Oxford. Moberg, A. 2013. Comparisons of simulated and observed Northern Hemisphere temperature variations during the past millennium – selected lessons learned and problems encountered: Tellus B 65, 19921, doi: 10.3402/tellusb.v65i0.19921. Nieto S. and Rodríguez-Puebla C. 2006. Comparison of Precipitation from observed data and general circulation models over the Iberian Peninsula. Journal of Climate, 19: 4254-4275. North, G. R., Cahalan, R. F., Coakley, J. A., 1981. Energy balance climate models. Review of Geophysica and Space Physics, 19:91-121, doi: 10.1029/RG019i001p00091 Otto, A., Otto, F., Boucher, O., Church, J., Hegerl, G., Forster, P. M., Gillett, N. P., Gregory, J., Johnson, G. C., Knutti, R., Lewis, N., Lohmann, U., Marotzke, J., Myhre, G., Shindell, D., Stevens, B., Allen, M.R. 2013. Energy budget constraints on climate response, Nature Geoscience, 6, 415– 416, doi:10.1038/ngeo1836 Peacock, J. A. 1983. Two-dimensional goodness-of-fit testing in astronomy, Monthly Notices of the Royal Astronomical Society, 202, 615-627, doi: 10.1093/mnras/202.3.615. Pelly, J. L., Hoskins, B. J. 2003. How well does the ECMWF Ensemble Prediction System predict blocking?. Quarterly Journal of the Royal Meteorological Society, 129, 1683–1702. doi: 10.1256/qj.01.173 Pennell, C., Reichler, T., 2011. On the Effective Number of Climate Models. Journal of Climate, 24, 2358–2367, doi: 10.1175/2010JCLI3814.1 Perkins, S. E., Pitman, A. J., N. H. Holbrook, McAneney, J., 2007. Evaluation of the AR4 climate models' simulated maximum temperature, minimum temperature and precipitation over Australia using probability density functions. Journal of Climate 20, 4356-4376, doi: 10.1175/JCLI4253.1 40 Radić, R., Clarke, G. K., 2011. Evaluation of IPCC Model's performance in simulating LateTwentieth-Century climatologies and weather patterns over North America. Journal of Climate 24, 5257-5274, doi: 10.1175/JCLI-D-11-00011.1 Randall, D. A., Wood, R. A., Bony, S., Colman, R. Fichefet, F., Fyfe, J., Kattsov, V., Pitman, A., Shukla, J., Srinivasan, J., Stouffer, R. J., Sumi, A., Taylor, K. E., 2007. Climate models and their evaluation. In: Solomon, S., Qin, D., Manning, Chen, Z., M., Marquis, M., Averyt, K., Tignor, M. M. B., Miller, H. L., Climate Change 2007, The Physical Science Basis, Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge , UK and New York, NY, USA, 996pp. Reichler, T. and Kim, J. K., 2008. How well do coupled models simulate today's climate? Bulletin of the American Meteorological Society 89, 303-311, doi: 10.1175/BAMS-89-3-303 Reifen, C. and Toumi, R., 2009. Climate projections: past performance no guarantee of future skill? Geophysical Research Letters, 36, L13704, doi: 10.1029/2009GL038082. Russell J., Stouffer R.J. and Dixon K.W. 2006. Intercomparison of the Southern Ocean circulations in IPCC coupled model control simulations. Journal of Climate, 19: 4560–4575. Saha, S., Moorthi, S., Pan, H., Wu, X., Wang, J., Nadiga, S., Tripp, P., Kistler, R., Woollen, J., Behringer, D., Liu, H., Stokes, D., Grumbine, R., Gayno, G., Wang, J., Hou, Y., Chuang, H., Juang, H. H., Sela, J., Iredell, M., Treadon, R., Kleist, D., Van Delst, P., Keyser, D., Derber, J., Ek, M., Meng, J., Wei, H., Yang, R., Lord, S., Van Den Dool, H., Kumar, A., Wang, W., Long, C., Chelliah, M., Xue, Y., Huang, B., Schemm, J., Ebisuzaki, W., Lin, R., Xie, P., Chen, M., Zhou, S., Higgins, W., Zou, C., Liu, Q., Chen, Y., Han, Y., Cucurull, L., Reynolds, R. W., Rutledge, G., Goldberg, M., 2010. The NCEP Climate Forecast System Reanalysis. Bulletion of the American Meteoological Society, 91, 1015–1057, doi: 10.1175/2010BAMS3001.1 Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., Behringer, D., Hou, Y., Chuang, H., Iredell, M., Ek, M., Meng, J., Yang, R., Mendez Malaquías, P., van den Dool, H., Zhang, Q., Wang, W., Chen, M., Becker, E., 2014. The NCEP Climate Forecast System Version 2. Journal of Climate, 27, 2185–2208. doi: 10.1175/JCLI-D-12-00823.1 41 Santer, B. D., Mears, C., Doutriaux, C., Caldwell, P., Gleckler, P. J., Wigley, T. M. L., Solomon, S., Gillet, N. P., Ivanova, D., Karl, T. R., Lanzante, J. R., Meehl, G. A., Stott, P. A., Taylor, K. E., Thorne, P. W., Wehner, M. F., Wentz, F. J., 2011. Separating signal and noise in atmospheric temperature changes: The importance of timescale. Journal of Geophysical research, 116, D22105, doi: 10.1029/2011JD016263. Schwalm, C. R., Huntinzger, D. N., Michalak, A. M., Fisher, J. B., Kimball, J. S., Mueller, B., Zhang, K. and Zhang, Y., 2013. Sensitivity of inferred climate model skill to evaluation decisions: a case study using CMIP5 evapotranspiration. Environmental Research Letters, 8, 024028, doi: 10.1088/1748-9326/8/2/024028 Scott, D. W. 1992. Multivariate Density Estimation: Theory, Practice and Visualization, John Wiley and Sons, New York, 336 pp, doi: 10.1002/9780470316849 Silverman, B. W., 1986. Density Estimation for Statistics and data Analysis, Chapman and Hall, London, 175pp. Stewart, D. E. and Leyk, Z., 1994. Meschach: Matrix Computations in C, Centre for Mathematics and its Applications, the Australian National University, Canberra, Australia. Stocker, T.F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M. (eds.). 2013. IPCC, 2013: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp. Sundberg, R., Moberg, A., and Hind, A. 2012. Statistical framework for evaluation of climate model simulations by use of climate proxy data from the last millennium – Part 1: Theory, Climate of the Past, 8, 1339-1353, doi: 10.5194/cp-8-1339-2012. Taylor, K. E., 2001. Summarizing multiple aspects of model performance in a single diagram, Journal of Geophysical Research, 106, 7183–7192, doi: 10.1029/2000JD900719. Taylor, K. E., Stouffer, R. J., Meehl, G. A., 2012. An overview of CMIP5 and the experiment