A new efficiency metric for the spatial evaluation and inter-comparison of climate and geoscientific model output

Developing and evaluating spatial efficiency metrics is essential for assessing how well climate or other models of the Earth's system reproduce the observed patterns of variables like precipitation, temperature, atmospheric pollutants, and other environmental data presented in a gridded format. In this study, we propose a new metric, the Modified Spatial Efficiency (MSPAEF), designed to overcome limitations identified in existing metrics, such as the Spatial Efficiency (SPAEF), the Wasserstein Spatial Efficiency (WSPAEF), or the Spatial Pattern Efficiency metric (E_sp). The performance of MSPAEF is systematically compared to these metrics across a range of synthetic data scenarios characterized by varying spatial correlation coefficients, biases, and standard deviation ratios. Results demonstrate that MSPAEF consistently provides robust and intuitive performance, accurately capturing spatial patterns under diverse conditions. Additionally, two realistic but synthetic case studies are presented to further evaluate the practical applicability of the metrics. In both examples, MSPAEF delivers results that align with intuitive expectations, while the other metrics exhibit limitations in identifying specific features in at least one case. Finally, as a real-world application, we rank global Coupled Model Intercomparison Project phase 6 (CMIP6) model data according to their skill in representing precipitation and temperature using the four different metrics. This application highlights that the MSPAEF rankings are most similar with E_sp with a normalized absolute ranking difference of 2.8 for precipitation, and 3.8 for temperature. These findings highlight the added value of the MSPAEF metric in evaluating spatial distributions and its potential to be used in climate or other environmental model evaluation or inter-comparison exercises.

1 Introduction

An accurate and comprehensive evaluation of climate models is crucial for understanding their limitations, enhancing the representation of processes, and distilling the most credible information on regional climate changes, impacts, and risks (Eyring et al., 2019). Climate models and their output are essential tools for scientists, practitioners and policymakers. However, before any application, it is important to use robust validation metrics to assess their performance and accuracy (Wagener et al., 2022). Traditional evaluation methods often emphasize temporal correlations or mean biases, while they often overlook crucial spatial patterns that significantly affect the model performance assessment. Capturing these spatial distributions is particularly important for environmental variables like precipitation and temperature, as the spatial variation greatly influences regional processes and subsequent applications, such as impact assessments (Citrini et al., 2024; Ramirez-Villegas et al., 2013).

In this paper, the term metric is used in a broad sense to refer to all indicators, statistics or distance measures that take as input two datasets and output a value that quantifies their relative performance or similarity. This includes, but is not limited to, quantities that satisfy the formal definition of a metric.

Many studies focus their evaluation of climate models on more simple metrics such as the Normalized Mean Square Error (NMSE) (Simpson et al., 2020), Bias (Mehran et al., 2014; Su et al., 2013), the Root Mean Square Error (RMSE) (Sillmann et al., 2013; Kamworapan and Surussavadee, 2019; Srivastava et al., 2020; Li et al., 2021; Nishant et al., 2022; Kim et al., 2020b), the Mean Absolute Error or Deviation (MAE or MAD) (Zittis and Hadjinicolaou, 2017; Lovino et al., 2018), the Relative Error (RE) (Kim et al., 2020a), or the Skill Score (Srivastava et al., 2020). For instance, for the evaluation and inter-comparison of downscaling approaches for climate models, Maraun et al. (2015) used a variety of metrics such as the Relative Error (RE) for the spatial evaluation of different downscaling methodologies, while for the temporal evaluation, performance indicators such as the mean square error (MSE), the correlation coefficient, the bias and Relative Error were used. A different concept was introduced by Brands et al. (2011) who evaluated global General Circulation Models (GCMs) using the overlap of the Probability Density function of the time series for each grid box between the model and the observational dataset. The Kolmogorov–Smirnov test (KS test) is another metric that was used for the evaluation of models contributing to the Coupled Model Inter-comparison Project (CMIP) (Brands et al., 2013). Kotlarski et al. (2014) evaluated a large ensemble of Regional Climate Models (RCM) contributing to the EURO-CORDEX initiative using a multitude of metrics, including spatial efficiency ones. This selection includes the spatial pattern correlation coefficient of time-averaged values (PACO), and the ratio of the spatial standard deviation (RSV) of time-averaged values between the model and the reference dataset. The Pearson Correlation Coefficient can also be used in the time series of seasonal means of a variable of each grid-box of the data, to assess how well a climate model represents the local location representativeness (Maraun and Widmann, 2015). While these metrics can be effective for certain applications, their performance can vary significantly based on the statistical and spatial characteristics of the data. For instance, the Pearson correlation coefficient is useful for capturing the linear relationship between two datasets; however, it does not account for systematic biases in the mean or differences in variance (i.e., scale differences). Similarly, other metrics like RMSE and MAE are sensitive to both magnitude and distribution of errors, but they may be disproportionately affected by a large bias in the mean, which can lead to a misrepresentation of spatial patterns. In contrast, the Kolmogorov–Smirnov test compares the underlying distributions of two datasets but lacks spatial context, which is often crucial in geoscientific modeling. These distinct limitations underscore the challenges of using traditional metrics on variables with uneven spatiotemporal distributions, like precipitation, which is often sparsely distributed in space. They also complicate comparisons between models that possess differing statistical properties, such as varying means or levels of variation.

Compound metrics that integrate several simpler approaches have been developed to assess the performance of geoscientific models more accurately. A commonly used goodness-of-fit indicator is the Kling–Gupta Efficiency (KGE), which is a robust compound metric, commonly used in hydrological sciences for comparing simulations to observations (Gupta et al., 2009; Deepthi and Sivakumar, 2022; Zittis et al., 2017). The SPAtial EFficiency metric (SPAEF) proposed by Demirel et al. (2018b, a) was inspired by KGE (Koch et al., 2018) and has emerged as a promising tool for evaluating hydrological models by considering three key aspects of spatial pattern accuracy: correlation coefficient, relative variability, and distribution of values in standardized space. This metric has also been applied to climate model evaluation, for the ranking of the performance of precipitation and temperature of different climate models (Lei et al., 2023; Deepthi and Sivakumar, 2022; Ahmed et al., 2019; Verma et al., 2023), as well as for the evaluation of precipitation in reanalyses datasets (Gomis-Cebolla et al., 2023). SPAEF combines various elements, equally weighted, into a single metric, offering a more comprehensive assessment of spatial distribution compared to traditional methods like mean-squared error or correlation alone. While SPAEF can be useful for evaluating climate models, it has some limitations, for example, it is insensitive to the magnitude of biases. Although this feature can be advantageous in certain contexts, an assessment of biases is crucial for a complete evaluation of a model and should be taken into account.

An existing modification and improvement of SPAEF is the Wasserstein SPAEF (WSPAEF) (Gómez et al., 2024). Unlike the original version, WSPAEF is sensitive to model biases. However, its sensitivity is based on the absolute magnitude of the bias, making it inherently dependent on the units of the variable used. This means that the metric's response to bias may vary across different variables or datasets with different scales, which should be considered when interpreting results. Another common limitation shared by both approaches is that two out of their three components – the standard deviation and differences in histogram values for SPAEF and the standard deviation and Wasserstein distance for WSPAEF – primarily describe the overall distribution of the parameter under evaluation rather than its spatial structure. This leaves only one component that evaluates the model variable based on the pattern of the variable, which is the spatial correlation coefficient.

Another indicator that was inspired by SPAEF is the Spatial Pattern Efficiency Metric (E_sp) (Dembélé et al., 2020b). This metric, like SPAEF, is bias-insensitive, and focuses on the spatial distribution of the patterns of a variable. This metric has been used to evaluate the spatial patterns of the hydrological processes in gridded precipitation datasets (Dembélé et al., 2020a). The SPAEF metric and some of its proposed modifications have been comprehensively evaluated by Yorulmaz et al. (2024).

Here, in response to the need for a comprehensive and multi-faceted evaluation of climate and other earth system models, we propose a new metric. This is based on a modification of SPAEF and is designed to more precisely capture the spatial distribution characteristics of climate model output. This new approach addresses certain limitations of existing metrics by incorporating spatial components that improve sensitivity to the spatial distribution and the relative bias of the variable under evaluation. We apply this new metric to synthetic datasets that imitate statistical properties and possible distributions of two types of climate variables, demonstrating its ability to offer improved insight into the spatial fidelity of model outputs. Through this work, we aim to provide the climate modeling community with a refined tool for assessing spatial distribution accuracy, contributing to more reliable model evaluations, improved climate model development, and, ultimately, more accurate projections. The proposed approach holds significant potential for applications in other scientific fields, including hydrology and environmental sciences.

It is important to stress that the MSPAEF metric is not intended to replace traditional multi-metric evaluation approaches. Instead, it is designed to complement existing measures by providing a balanced indicator that captures both spatial pattern similarity and relative mean bias within a single metric. In this way, MSPAEF serves as a useful additional tool, particularly when a unified summary statistic is desirable.

2 Data and Methods

2.1 Metrics

This study focuses on compound metrics specifically designed for the spatial evaluation of climate or other environmental parameters. The four metrics evaluated and inter-compared in this work are presented below.

2.1.1 SPAEF metric

The SPAEF (SPAtial EFficiency) metric by Demirel et al. (2018b) is a robust metric that was originally created for hydrological model evaluation, and it is used to characterize the performance of a model regarding the spatial distribution of a variable. It is defined as:

$$\begin{array}{}\text{(1)}& \text{SPAEF}=\mathrm{1}-\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\gamma }-\mathrm{1}{)}^{\mathrm{2}}}\end{array}$$

with:

$$\begin{array}{}\text{(2)}& \mathit{\alpha }=\mathit{\rho }(M,O),\mathit{\beta }={\displaystyle \frac{\left(\frac{{\mathit{\sigma }}_{\mathrm{M}}}{{\mathit{\mu }}_{\mathrm{M}}}\right)}{\left(\frac{{\mathit{\sigma }}_{\mathrm{O}}}{{\mathit{\mu }}_{\mathrm{O}}}\right)}}\text{ and }\mathit{\gamma }={\displaystyle \frac{{\sum }_{j=\mathrm{1}}^n\text{min}(K_j,L_j)}{{\sum }_{j=\mathrm{1}}^n{K}_j}}\end{array}$$

where ρ is the spatial correlation coefficient between the points of the model and the observations, σ and μ are the standard deviation and mean value for each of the model and observation, K and L are histograms with n common bins, of the standardized values (z-score) of the model and observations, respectively, and γ is the histogram overlap of the standardized values of the two datasets. This metric is bias-insensitive since none of its three terms are affected by present bias, and therefore the final metric conveys only the similarity in the patterns between the model and the observations. The SPAEF metric takes values from −∞ to 1, with 1 indicating an excellent agreement of the model with the observations.

2.1.2 The Spatial Pattern Efficiency Metric (E_sp)

The E_sp metric by Dembélé et al. (2020b) was designed to be a simple and robust metric for hydrological model evaluation, and it is defined as:

$$\begin{array}{}\text{(3)}& E_{\text{sp}}=\mathrm{1}-\sqrt{(r_{\mathrm{s}}-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\gamma }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}}\end{array}$$

with:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{r}_{\mathrm{s}}& =\mathrm{1}-{\displaystyle \frac{\mathrm{6}{\sum }_{\mathrm{1}}^n{d}^{\mathrm{2}}}{n(n^{\mathrm{2}}-\mathrm{1})}},\mathit{\gamma }={\displaystyle \frac{\left(\frac{{\mathit{\sigma }}_{\mathrm{M}}}{{\mathit{\mu }}_{\mathrm{M}}}\right)}{\left(\frac{{\mathit{\sigma }}_{\mathrm{O}}}{{\mathit{\mu }}_{\mathrm{O}}}\right)}}\text{ and }\\ \mathit{\alpha }& =\mathrm{1}-E_{\text{RMS}}(Z_{\mathrm{M}},Z_{\mathrm{O}})\end{array}\end{array}$$

where r_s is the Spearman rank-order correlation coefficient, σ and μ are the standard deviation and mean value for each of the model and observation, Z_M and Z_O are the standardized values (z-score) of the model and observations, respectively, and E_RMS is the root mean square error of the standardized values. A benefit of this metric is that it does not require any user-defined parameters, such as the bins of the histograms, and that two of the three components take into account the spatial distribution of the variable. It takes values between −∞ and 1, with 1 indicating excellent performance of the model.

2.1.3 WSPAEF Metric

An relatively newly developed modification of the SPAEF is the Wasserstein Spatial Efficiency (WSPAEF) by Gómez et al. (2024) and is defined as:

$$\begin{array}{}\text{(5)}& \text{WSPAEF}=\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\sigma }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\varphi }{)}^{\mathrm{2}}}\end{array}$$

with:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\alpha }& =\mathit{\rho }(M,O),\mathit{\sigma }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{M}}}{{\mathit{\sigma }}_{\mathrm{O}}}}\text{ and }\\ \mathit{\varphi }& =\text{WD}={\left(\sum _{i=\mathrm{1}}^n|X_{\left(i\right)}-Y_{\left(i\right)}{|}^p\right)}^{\frac{\mathrm{1}}{p}}\end{array}\end{array}$$

where ρ is the spatial correlation coefficient between the points of the model and the observations, σ is the standard deviation for each of the model and observation, WD is the Wasserstein distance of order p=2, with X_(i) and Y_(i) being the ith order statistic of the samples of the observations and model, respectively, and n is the total number of samples in each dataset. This means that the values of the two datasets have been arranged in ascending order, and X_(i) is the ith value in the sorted list. WD was calculated using the original values of the two datasets to explicitly account for the bias. This modification of SPAEF is sensitive to bias since the Wasserstein distance term is a measure of the minimum required effort that is needed to change the histogram of one dataset to match the histogram of the other dataset. It takes values from 0 to ∞, with values approaching zero indicating excellent agreement of the model with the observations.

2.1.4 MSPAEF Metric

According to these definitions, some limitations can be introduced when using the three previous metrics. For instance, SPAEF and WSPAEF only include one component that considers the spatial distribution of the variable's values, which is the correlation coefficient in the space domain, while their other two components take into consideration the overall distribution of the parameter. Furthermore, SPAEF and E_sp are bias-insensitive metrics. While this characteristic may be beneficial in certain contexts, it is not ideal for a comprehensive assessment of climate models, where capturing both spatial patterns and biases is essential for a thorough evaluation. In addition, all three of them are scale-dependent, meaning that their values change based on the units of the input variable.

Therefore, we propose a new metric, on the basis of SPAEF that is bias-sensitive, but at the same time scale-independent. This metric conveys information not only about the spatial patterns but also the biases, like the WSPAEF. This is desired because, for a more complete evaluation of the spatial characteristics of a modeled parameter, both the pattern information, and the bias are essential. This metric, like SPAEF and E_sp, takes values from negative infinity to 1, with 1 indicating a perfect agreement of the model compared to the observations. The proposed modification to the SPAEF is the following:

$$\begin{array}{}\text{(7)}& \text{MSPAEF}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{4}}}}\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }{)}^{\mathrm{2}}+(\mathit{\gamma }{)}^{\mathrm{2}}+(\mathit{\delta }{)}^{\mathrm{2}}}\end{array}$$

with:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{\alpha }& =\mathit{\rho }(\mathrm{M},\mathrm{O}),\mathit{\beta }=\text{NRMSE},\mathit{\gamma }={\displaystyle \frac{|\stackrel{\mathrm{‾}}{M}-\stackrel{\mathrm{‾}}{O}|}{{\text{IQR}}_{\mathrm{o}}}},\\ \mathit{\delta }& =\sqrt{{\left({\mathit{\sigma }}^{\mathrm{2}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}-\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}}}\right)}^{\mathrm{2}}}\end{array}\end{array}$$

also:

$$\begin{array}{}\text{(9)}& \text{NRMSE}={\displaystyle \frac{\mathrm{1}}{{\text{IQR}}_{\mathrm{O}}}}\sqrt{{\displaystyle \frac{{\sum }_i^n(M_i-O_i{)}^{\mathrm{2}}}{n}}}\text{ and }\mathit{\sigma }={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{M}}}{{\mathit{\sigma }}_{\mathrm{O}}}}\end{array}$$

where ρ is the spatial correlation coefficient between the grid cells of the model and the observations, NRMSE is the Root Mean Square Error calculated with the original data and normalized with the inter-quartile range of the observational data (IQR_o), while σ_O and σ_M are the standard deviations of the observations, and the model, respectively.

This modification of the SPAEF metric is sensitive to the relative value of bias, while at the same time capturing the spatial pattern differences. Specifically, the correlation coefficient (α term) ensures that the metric reflects the spatial agreement between model output and observations. The normalized RMSE (β term) is computed using the interquartile range (IQR) of the observations rather than the standard deviation, as IQR provides a more robust measure of variability (Rousseeuw and Hubert, 2011; Huber and Ronchetti, 2009). This term captures both spatial pattern agreement and bias, behaving differently depending on the magnitude of bias. When bias is small, it emphasizes spatial differences, whereas, in cases of large bias, it is primarily driven by magnitude differences. The relative bias term (γ term) is also normalized by the interquartile range of the observations and is introduced to explicitly quantify systematic differences in the mean. Finally, the δ term accounts for differences in the spread of values by incorporating the standard deviation ratio between the model and observations. This final term comprises two components to ensure symmetry around a value of 1, thereby ensuring the same result regardless of whether the standard deviation ratio or its reciprocal is applied.

An advantage of MSPAEF over SPAEF is that its definition does not rely on user-defined parameters. The MSPAEF metric is formulated without such arbitrary choices, unlike SPAEF, which requires user-specified bins for the histogram overlap term. This makes it more objective and consistent across different datasets, reducing the influence of subjective parameter selection on the results.

A further advantage of normalizing all components by IQR and expressing each term in dimensionless form is that MSPAEF becomes relatively insensitive to preprocessing decisions such as unit conversions, domain rescaling, or masking of small regions. Unlike WSPAEF, whose values are affected by the unit of choice, MSPAEF can produce comparable values under a wide range of analysis settings, which increases reproducibility across studies and domains.

An important motivation for this definition of MSPAEF is to penalize inconsistent model behavior, where a model with poor spatial correlation but a very small bias (or vice versa) can appear artificially good in composite metrics. Since MSPAEF treats each discrepancy in either the spatial patterns or in the mean relative bias as an orthogonal dimension, it prevents a strong performance in one characteristic from masking deficiencies in the other. This provides a more complete evaluation of how each model differs from observations.

Based on its definition, MSPAEF is sensitive to both spatial pattern agreement and mean bias, and this sensitivity can be directly interpreted in the corresponding spatial maps. Each component of the metric can be associated with an observable feature in the maps. The spatial correlation term represents how well the locations and intensities of spatial gradients align between the model and observations. The NRMSE identifies differences in grid-point magnitudes across the domain that are linked to both the location and intensity of gradients, but also to the mean bias. The relative mean bias term quantifies systematic offsets in the average values of the field. Finally, the variability term reflects how well the model reproduces the amplitude of spatial fluctuations. Together, these components allow changes in the MSPAEF value to be directly linked with spatial mismatches.

2.2 Synthetic Data generation

For the development of MSPAEF and the comparison with existing metrics we used a methodology inspired by Gómez et al. (2024). This approach uses synthetic data for both the observations and the model, in order to identify the behavior of the metrics for a variety of predefined combinations of correlation, bias and standard deviation ratio between the two sets of data.

A 10 × 10 matrix of random data that incorporates some spatial correlation between the different grid points was created using a covariance model. To achieve that, a distance matrix was first created for the specific grid size, which holds the Euclidean distance of all possible pairs of points. Then, the covariance matrix was calculated using this distance matrix and the Mattern covariance function:

$$\begin{array}{}\text{(10)}& M_{\mathrm{v}}\left(\mathbf{r}\right)={\displaystyle \frac{{\mathrm{2}}^{\mathrm{1}-v}}{\mathrm{\Gamma }\left(v\right)}}\times {\displaystyle \frac{{\mathbf{r}}_{\mathrm{s}}^v}{l_{\mathrm{s}}}}\times K_{\mathrm{v}}\left({\mathbf{r}}_{\mathrm{s}}\right)\end{array}$$

with:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{l}_{\mathrm{s}}& =\mathrm{1.0},v=\mathrm{1.5},\mathrm{\Gamma }\left(v\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{t}^{v-\mathrm{1}}{e}^{-t}\mathrm{d}t\text{ and }\\ {\mathbf{r}}_{\mathrm{s}}& =\sqrt{\mathrm{2}\times v\times \mathbf{r}}\end{array}\end{array}$$

where v is the smoothing parameter, r is the distance matrix, r_s is the scaled distance matrix, l_s is a scaling factor, Γ(v) is the gamma function and K_v is the modified Bessel function of the second kind.

To create the synthetic field with prescribed properties such as a specific standard deviation and mean value, the multivariate normal function was used. This function requires two inputs: a mean vector μ that specifies the expected value for each point on the grid and a covariance matrix Σ encoding the relationships (e.g., variances and covariances) between the values at all pairs of points on the grid. For the mean vector, a vector that holds the target mean values was created. As a covariance matrix input, the following matrix was used in order to achieve the desired standard deviation:

$$\begin{array}{}\text{(12)}& C={\left({\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{M}}}}\right)}^{\mathrm{2}}\times M_{\mathrm{v}}\left(\mathbf{r}\right)\end{array}$$

with

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_{\mathrm{M}}=\sqrt{{\displaystyle \frac{\text{Tr}\left(M_{\mathrm{v}}\right(\mathbf{r}\left)\right)}{d}}}\end{array}$$

where σ is the target standard deviation, M_v(r) is the Mattern Covariance function calculated above and d is the number of spatial points. The output of the multivariate normal function represents the spatial distribution of the observational dataset.

2.3 Correlated data generation

A second gridded dataset was then created with some correlation, standard deviation ratio and bias target compared to the first grid, to simulate the spatial distribution of the model dataset, using the following equation:

$$\begin{array}{}\text{(14)}& y={\mathit{\lambda }}_{\mathrm{t}}\times \left({\mathit{\rho }}_{\mathrm{t}}\times (\mathbf{x}-\stackrel{\mathrm{‾}}{x})+\sqrt{\mathrm{1}-{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{2}}}\times \mathbf{A}\times {\mathit{\sigma }}_x\right)+\stackrel{\mathrm{‾}}{x}+{\mathit{\delta }}_{\mathrm{t}}\end{array}$$

where λ_t is the standard deviation ratio target, ρ_t is the correlation coefficient target, δ_t is the bias target, x is the original matrix, $\stackrel{\mathrm{‾}}{x}$ is the mean value of the original matrix x, σ_x is the standard deviation of the original matrix x and A is a matrix of the same dimensions as x, filled with random numbers drawn from the normal distribution.

2.4 Skewed data generation

The above data are meant to represent climate variables that have a normal distribution since the random data used follow a normal distribution. Nevertheless, numerous climate and other geoscientific model output variables do not follow a normal distribution, but instead exhibit skewed or exponential distributions, as is the case with daily precipitation (Ensor and Robeson, 2008). To simulate skewed distributions, the normally distributed observational values were transformed using exponentiation. A corresponding model dataset was generated by first establishing a specific correlation and standard deviation ratio between the underlying normally distributed values, and then introducing the desired bias after the exponentiation transformation was applied to both datasets. A limitation of this approach is that even though the target correlation coefficient is generally well-preserved through the exponentiation, the target standard deviation ratio and bias are not directly translated to the skewed data. Nevertheless, relative changes in the target standard deviation ratio and bias are reflected in the resulting skewed data.

Using the original observational grid defined above, the skewed observations and model grids are defined as:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{x}_{\mathrm{s}}& =e^{{\mathbf{x}}_{\text{rm}}}\\ y_{\mathrm{s}}& =e^{{\mathit{\lambda }}_{\mathrm{t}}\times \left({\mathit{\rho }}_{\mathrm{t}}\times \left({\mathbf{x}}_{\text{rm}}\right)+\sqrt{\mathrm{1}-{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{2}}}\times \mathbf{A}\times {\mathit{\sigma }}_x\right)}+{\mathit{\delta }}_{\mathrm{t}}\end{array}\end{array}$$

with:

$$\begin{array}{}\text{(16)}& {\mathbf{x}}_{\text{rm}}=\mathbf{x}-\stackrel{\mathrm{‾}}{x}\end{array}$$

where x is the original matrix, x_rm is the original matrix after removing the mean, and the other terms are the same as defined in Eq. (14).

Figure 1 Examples of synthetic data generation. Panels (a) and (b) show the spatial distribution and probability density function (PDF), respectively, of synthetic data generated with a normal distribution. Panels (c) and (d) show the corresponding spatial distribution and PDF for skewed data.

Although the data created this way follow a highly positively skewed distribution, their shape closely resembles an exponential distribution when visualized with an insufficient number of histogram bins, due to the generation process, as seen in Fig. 1.

The four metrics were calculated for many combinations of bias, standard deviation ratio and correlation, for both the normal and skewed distribution cases. For each combination of the aforementioned parameters, and for both the normally and skewed distributed cases, the procedure was repeated 200 times, and the median value of each metric was used to ensure convergence of the metrics. While the synthetic data are generated to match specific target values of correlation, bias, and standard deviation ratio between them, the stochastic nature of the process means these targets are only approximately achieved. Additionally, non-trivial spatial variations may still occur across realizations, which can affect different metrics in distinct ways. The repetition and the use of the median metric values help reduce the influence of these variations and provide a more robust estimate of each metric's behavior under the intended conditions. To facilitate an easier comparison of the performance of the metrics, these were modified slightly so that they take values equal to or greater than zero, with zero indicating perfect agreement with the observations (see Appendix A).

Figure 2 Examples of metric behaviour curves for a range of correlation and bias values. Black and purple curves indicate a well-behaved metric, while the blue curves indicate a poorly behaving metric. In the figure, λ is the standard deviation ratio target, δ is the bias target, and ρ is the correlation target between the model and the observations.

In these adjusted conditions, where zero indicates a perfect match with the observations, a well-behaved metric is expected to show decreasing values as the correlation coefficient increases and the bias decreases. In the λ-δ-ρ plots, this is reflected as the curves shift at lower coordinates values as one moves towards the right and upward parts of the subplots, as illustrated by the purple curves in Fig. 2. Additionally, the lowest metric value for any combination of correlation and bias is anticipated at a standard deviation ratio of 1. This is reflected in the curve minimum being at a standard deviation ratio of 1 for each subplot, with metric values increasing as the standard deviation ratio deviates from 1.

In the examples of the top left subplot of Fig. 2, the black and purple curves indicate well-behaving metrics, since for all of them, the minimum values are found at λ of 1, and the values increase monotonically as the standard deviation ratio deviates from 1, even though some of the curves are not perfectly symmetric. The blue curves indicate poorly behaving metrics because, for curve e, the minimum value is not found at a standard deviation ratio of 1, while for curve d, the curve does not monotonically increase to the left side of the minimum.

2.5 Real climate data

The sixth phase of the Coupled Model Intercomparison Project (CMIP6) is the latest generation of CMIP, which is designed to improve understanding of past, present, and future climate change driven by greenhouse gas emissions. Numerous research institutions worldwide contributed their climate models, following standardized protocols for radiative forcing and data formatting. CMIP6 models have generally shown improved performance in simulating various climate fields compared to previous CMIP generations (Li et al., 2020). They also exhibit reduced biases in tropical regions, such as the double ITCZ bias (Tian and Dong, 2020). In this study, 33 CMIP6 models were evaluated against the ERA5 reanalysis dataset (Hersbach et al., 2020). For each model, the ensemble mean of all available variants was used, for precipitation and 2 m temperature during the historical period from 1981 to 2010. The CMIP6 models and the ERA5 dataset were remapped to 1° resolution using linear interpolation for the temperature variable and first-order conservative interpolation for the precipitation variable. For precipitation, the mean annual total in millimeters was calculated over the period, while for temperature, the mean annual value in Kelvin was computed.

To compute MSPAEF or any other spatial metrics, the model and observational datasets need to be placed into a common grid. Regridding data can generally affect the values of metrics that rely on grid-point values. However, the relative performance rankings among models are expected to remain largely unchanged because all models undergo the same regridding procedure to the same target grid. Therefore, any smoothing or distortion that occurs due to the interpolation is expected to have a similar and small effect across the datasets, provided that the regridding does not involve extreme changes to resolution.

3 Results and Discussion

In Sects. 3.1 and 3.2, the MSPAEF, SPAEF and E_sp metrics were modified slightly so that their best performance is indicated at the zero value, like in the WSPAEF metric, for easier inter-comparison of their performance. Contrariwise, in Sect. 3.3, all metrics were applied in their original form.

3.1 Metrics Behavior

For each of the four metrics examined, we calculated its value for a range of values of spatial correlation coefficient, bias and standard deviation ratio between the observational data and the model. The correlation coefficient, bias, and standard deviation ratio were each sampled at discrete intervals within the following ranges: correlation (−0.9 to 0.9), bias (0 to 3), and standard deviation ratio (0.3 to 1.8). The following plots show the behavior of the four metrics following the logic presented in Fig. 2, using a mean of 10 and a standard deviation of 1 for the observational data.

Figure 3 SPAEF metric behavior (λ is the standard deviation ratio target, ρ is the correlation target and δ is the bias target between the model and the observation dataset).

For the SPAEF metric (Fig. 3), the values generally decrease with increasing correlation, which is the expected behavior of a well-defined spatial efficiency metric. In addition, a good behavior is observed for each curve in the cases where there is no bias present, as the minimum of the curves in these cases is found at standard deviation ratios of 1. Nevertheless, as the bias increases ($\mathit{\delta }=\mathrm{1.0},\mathrm{3.0}$), the minimum of the curves generally shifts to values of a standard deviation ratio greater than 1. Specifically, for the skewed distributed variable, it moves for both cases where there is a bias present, while for the normally distributed variable, the shift of the minimum to standard deviation ratios values greater than 1 is observed only for the cases of a large bias (δ=3.0), with the shift of the normally distributed cases being more perceivable with a high correlation coefficient (ρ=0.9). This indicates that the metric does not perform as intended in larger biases, especially in the case of a high correlation coefficient, where this shift of the minimum is more evident. In addition, the increasing bias does not directly affect the magnitude of the SPAEF values, highlighting that this metric is bias-insensitive.

Figure 4 E_sp metric behavior (λ is the standard deviation ratio target, ρ is the correlation target and δ is the bias target between the model and the observation dataset).

The behavior of the Spatial Pattern Efficiency Metric (E_sp) (Fig. 4) is very similar to that of SPAEF. There is good performance in the case of no bias, across the values of the correlation coefficient, but it performs worse as the bias increases, especially for the case of a high correlation coefficient (ρ=0.9), as the minimum of the curves shifts to standard deviation values greater than 1. Specifically, for the cases of a small bias (δ=1.0) the metric performs well for the normally distributed data, but for the skewed data, the minimum of the curves shifts to standard deviations greater than 1. Conversely, for the cases of a large bias (δ=3.0), while the skewed variables have their minimum shifted to standard deviation ratios larger than 1 for all correlation coefficients, the normally distributed cases also have a shift of the minimum to standard deviations greater than one that is perceivable for the high correlation coefficient only (ρ=0.9). This metric, like in the case of SPAEF, is bias-insensitive, since different bias values do not directly reflect to changes in the values of the metric.

Figure 5 WSPAEF metric behavior (λ is the standard deviation ratio target, ρ is the correlation target and δ is the bias target between the model and the observation dataset).

On the other hand, the WSPAEF metric (Fig. 5) exhibits a distinct behavior compared to the aforementioned metrics. Similarly to the two previous metrics, a good performance is observed for both the normal and skewed distributed data in the absence of bias. However, for the normally distributed cases across all combinations of correlation coefficient and bias, the minimum value of the curve is consistently at a standard deviation value of 1. For the skewed data, in the cases where there is a small bias present (δ=1.0) and the correlation coefficient is near 0 ($\mathit{\rho }=-\mathrm{0.3},\mathrm{0.3}$), the minimum value of the curve shifts at standard deviations ratio slightly below 1. When a substantial bias is present (δ=3.0), the minimum value of the curve for the skewed distributed data shifts to a standard deviation ratio significantly below 1. Additionally, unlike SPAEF and E_sp, this metric demonstrates sensitivity to bias, as the metric values generally increase with increasing values of bias.

Figure 6 MSPAEF metric behavior (λ is the standard deviation ratio target, ρ is the correlation target and δ is the bias target between the model and the observation dataset).

The MSPAEF metric, illustrated in Fig. 6, demonstrates a robust performance across all combinations of correlation coefficients and biases for both the normal and skewed distributed variables. Nonetheless, under conditions of large bias (δ=3.0) and negative correlation coefficients (ρ = −0.9, −0.3), the response curve for normally distributed data exhibits a near-plateau around the standard deviation ratio of 1. This suggests a diminished sensitivity to variations in the standard deviation ratio within this specific range. Similar to WSPAEF, MSPAEF is sensitive to bias, with values generally increasing with increasing bias. This trend is particularly pronounced for cases with a large positive correlation coefficient (ρ=0.9). Generally, the metric values decrease with increasing correlation and decreasing bias, which aligns with the expected behavior of a well-behaved metric.

Figure 7 (a) Example of spatial pattern of synthetic observations and Models A and B in Example 1 (Observations were to represent a scenario of the mean annual precipitation in mm of a specific region. Model A represents a model with high correlation and a small relative bias compared to the observations, while Model B represents a model with negative correlation and no bias compared to the observations), (b) boxplots of the metric values for Models A and B. (SPAEF: Spatial Efficiency, WSPAEF: Wasserstein Spatial Efficiency, MSPAEF: Modified Spatial Efficiency, E_sp: Spatial Pattern Efficiency)

3.2 Synthetic data examples

Two examples are presented to illustrate the differences in the values and the interpretation of the four metrics, using normally distributed synthetic data as described in the Methods section. These synthetic datasets are designed to realistically resemble observations and climate model output for evaluation. The first example shown in Fig. 7a, presents the synthetic observations and two different models with predefined statistical properties. The observations were generated with a mean value of 200 and a standard deviation of 70, representing a scenario such as the mean annual precipitation in millimeters (mm) in a specific region of the planet. Although daily precipitation often follows a highly skewed or exponential distribution, annual averages can be closely approximated to a normal distribution, especially outside of polar regions (see Appendix B). Model A was created using a spatial correlation coefficient equal to 0.85, a standard deviation ratio of 1.1 and a domain-average bias of −10. In contrast, Model B was created using a correlation coefficient of −0.4, a standard deviation ratio of 1.0, and a domain-average bias of 0. While the use of negative spatial correlations might look unusual, especially at global scales, they can occur regionally, especially in precipitation fields due to known large biases near the equator, such as the double ITCZ bias (Ma et al., 2023).

Model A shows a spatial pattern that resembles the observations closely, with only minor differences. In contrast, Model B displays a significantly different spatial pattern for the variable under evaluation. Intuitively, Model A should be considered the better model due to its close alignment with the observations and its minimal relative bias. On the other hand, Model B should be regarded as the less skillful model because its spatial pattern diverges considerably from the observations, even though it does not exhibit significant bias when looking at the domain average.

In Fig. 7b, the boxplots show the distribution of the values of the four metrics from 20 different realizations of the synthetic data generated with the aforementioned parameters. The WSPAEF metric indicates that Model B has better performance compared to Model A, since it has the lower values for the metric, while the other three metrics conclude that Model A illustrates a model with better performance compared to Model B, as our intuition suggests. The bad performance of WSPAEF in this example is attributed to its very high sensitivity to the actual values of the bias due to the Wasserstein distance component, which overwhelms the contribution of the other two components on the final value of the metric.

Figure 8 (a) Example of the spatial pattern of observations and Models A and B in Example 2 (Observations were to represent a scenario of the annual temperature in degrees Celsius of a specific region. Model A represents a model with medium correlation and no bias compared to the observations, while Model B represents a model with high correlation and a large bias compared to the observations), (b) boxplots of the metric values for Models A and B. (SPAEF: Spatial Efficiency, WSPAEF: Wasserstein Spatial Efficiency, MSPAEF: Modified Spatial Efficiency, E_sp: Spatial Pattern Efficiency)

The second example is shown in Fig. 8. The observations of this example were generated with a mean of 6 and a standard deviation of 1, representing a scenario such as the annual temperature in a region of the planet. Model A was created using a spatial correlation coefficient of 0.6, a standard deviation ratio of 1.1 and a domain-average bias of 0. Model B was created using a spatial correlation coefficient of 0.8, a standard deviation ratio of 2.0 and a domain-average bias of 7.5. While this large bias example might seem unrealistic at global scales, it can occur at local scales and/or momentarily (McSweeney et al., 2015).

Model A spatial pattern exhibits some correlation compared to the original data but with no significant bias present. In contrast, Model B has a very similar spatial pattern for the variable, but with a very pronounced positive bias. Based on these observations, intuitively, Model A should be considered the better of the two models. Although its correlation with the observations is average, the lack of bias makes it more reliable. On the other hand, Model B should be considered the worst of the two because, despite its similarity to the observed spatial pattern, it has a substantial bias compared to the observations.

In Fig. 8b the boxplots show the distribution of the values of the four metrics from 20 different realizations of the synthetic data generated with the aforementioned parameters. The WSPAEF and MSPAEF metrics indicate that the Model A has far better performance compared to the Model B. This is expected, since there is a large bias present in Model B which dramatically increases the values of these two bias-sensitive metrics, contrariwise to Model A which has no bias present. On the other hand, the other two metrics indicate that Model B illustrates a model with better performance compared to Model A, because these two metrics are bias-insensitive and therefore are affected more by the underlying patterns of the variable in the two models.

Figure 9 Scatter plot of the Modified Spatial Efficiency (MSPAEF) values for annual mean precipitation and 2 m temperature derived from 33 CMIP6 models during the historical period of 1981–2010.

3.3 Application to real data

In this section, the four metrics were used in their original form, where the value of 1 indicates the best performance for MSPAEF, SPAEF and E_sp, while for WSPAEF the best performance is indicated at a value of 0. The metrics were applied to global CMIP6 model output, in order to benchmark them against the ERA5 reanalysis for annual precipitation totals and mean near-surface temperature. For both examples, we averaged over the 1981–2010 period. These multi-decadal averages were used to reduce the short term variability, and highlight the long term climatological signal (Du et al., 2022; Nooni et al., 2023). The performance of the CMIP6 models using the MSPAEF metric for the two variables is presented in Fig. 9. Some models perform very well for the spatial distribution of both variables (e.g., CESM2), while others perform well for one of the two variables, and much poorer for the other. For instance, the MIROC6 model is found to perform better than most other CMIP6 models for precipitation, but it is the second worst-performing model in representing the spatial distribution of near-surface temperature. On the other hand, ACCESS-ESM1-5 performs very well for the temperature, being in the top few models, but it is the worst-performing model in representing precipitation. Then, there are models that perform badly for both variables, such as MIROC-ES2L which performs below average in terms of the precipitation variable and is the worst-performing model by far in the representation of temperature.

In line with a previous study (Jun et al., 2008), we generally expect models from the same developers or model family to cluster together in performance, and we do observe this behavior for most groups. However, a few clear outliers emerge. These outliers typically correspond to models with substantially coarser atmospheric resolution (see Appendix C), which can markedly alter precipitation characteristics and therefore degrade overall performance relative to their higher-resolution counterparts. Although differences in parameterizations and internal variability may also play a role, the systematic link between resolution and performance offers the most plausible explanation for the departures from the expected clustering.

Regarding the actual values of the metric, all CMIP6 models except MIROC-ES2L, have values greater than 0.9 for the temperature variable. MRI-ESM2-0 is the model with the best performance for the temperature variable, which indicates that this model better represents not only the spatial patterns, but also the actual values of temperature at each grid-point. Contrariwise, all models except the two CESM2s and the NorESM2-MM models have values less than 0.85 for the precipitation variable. This indicates that the CMIP6 models generally can capture well the spatial distribution and magnitude of the temperature, while they struggle a lot more in the representation of the precipitation features.

Although comparing the MSPAEF values for two different variables is not straightforward, some conclusions on the nature of the variables can be deduced. The spatial distribution of temperature is influenced by factors such as orography, latitude, and proximity to the oceans, which are easier to be represented in climate models, due to their generally smoother and slow-varying nature in space. On the other hand, precipitation tends to have steeper gradients. Modeling precipitation is significantly more complex (Legates, 2014; Räisänen, 2007), as it is largely influenced by sub-grid-scale processes like convection and cloud microphysics that are not explicitly resolved in the models, and they need to be parametrized (Pieri et al., 2015), which makes it more difficult for climate models to correctly capture the underlying gradients and the actual grid-point values.

Table 1 Values of the four metrics and CMIP6 model rank for the annual total precipitation, averaged over the 1981–2010 period (SPAEF: Spatial Efficiency, WSPAEF: Wasserstein Spatial Efficiency, MSPAEF: Modified Spatial Efficiency, E_sp: Spatial Pattern Efficiency).

Table 2 Values of the four metrics and CMIP6 model rank for the monthly 2m temperature variable averaged over the 1981–2010 period (SPAEF: Spatial Efficiency, WSPAEF: Wasserstein Spatial Efficiency, MSPAEF: Modified Spatial Efficiency, E_sp: Spatial Pattern Efficiency).

Tables 1 and 2 present the values and the corresponding ranking of the models, based on the four metrics discussed in this study, for annual precipitation and 2 m temperature, respectively. The metrics generally exhibit a degree of consensus regarding the identification of the best and worst-performing models for each variable. However, the precise ranking order of the models varies across the different metrics.

For precipitation (Table 1), the model rankings derived from the WSPAEF metric differ the most from those obtained from the other three metrics, as illustrated in Fig. 10. The normalized sum of absolute rank differences between the CMIP6 models exceeds the 7 units when comparing the WSPAEF metric with the other three metrics. This significant differentiation in the WSPAEF rankings can be attributed to the very large values this metric attains as a consequence of its high sensitivity to the inherent biases that arise from its Wasserstein distance component. Given that annual precipitation was the variable under consideration, a bias of just 30 mm in the globally averaged total annual precipitation can result in WSPAEF metric values exceeding the 30 units, even if this bias represents only approximately 3 % of the total annual global precipitation.

Figure 10 Normalized sum of absolute rank differences between pairs of metrics for the 33 CMIP6 models, based on mean annual precipitation and temperature during the historical period of 1981–2010. (SPAEF: Spatial Efficiency, WSPAEF: Wasserstein Spatial Efficiency, MSPAEF: Modified Spatial Efficiency, E_sp: Spatial Pattern Efficiency).

The MSPAEF metric rankings exhibit the most similarity to those of E_sp with a rank difference value of 2.8 and the least similarity to WSPAEF with a rank difference value of 7.9, with SPAEF demonstrating an intermediary level of similarity. This suggests that MSPAEF effectively captures the spatial distribution of precipitation while at the same time not being overwhelmed by existing bias.

For temperature (Table 2), the four metrics generally exhibit a greater degree of agreement regarding the ranking of the models. However, more variability is observed in the identification of the top-performing models. Concerning the WSPAEF metric in Table 2, its values exceed the 2 units only for three models. This is attributed to the considerably smaller absolute values of the mean bias between the individual models and the observations, measured in Kelvin, which rarely exceed one or two degrees. This relatively small absolute value of the mean bias allows the other components of the WSPAEF metric, like the correlation coefficient, to exert a greater influence in the final metric value. This enables WSPAEF to achieve greater consistency in its rankings with the other three metrics for temperature than for precipitation, as depicted in Fig. 10.

The most substantial differences in temperature rankings are observed between WSPAEF and SPAEF, while the smallest difference occurs between the E_sp and SPAEF metric. The rankings of MSPAEF exhibit the greatest similarity to the E_sp metric with rank differences value of 3.8, and the least similarity to SPAEF with rank differences value of 5.5.

The similarity of MSPAEF to the bias-insensitive metrics might look unusual at first glance (compared to the other bias-sensitive metrics), given that MSPAEF is a bias-sensitive metric. However, this is less surprising when considering the large sensitivity of WSPAEF (the other bias-sensitive metric) to the absolute value of the mean bias. This sensitivity is what causes WSPAEF's performance to diverge significantly from that of MSPAEF (and the other two bias-insensitive metrics), particularly for the precipitation variable, where the bias in the mean can be large in absolute terms.

Using weights for the different components could serve as a way to improve the performance of the existing metrics, to closely match the performance of MSPAEF. For example, in the case of WSPAEF, achieving a behavior more consistent with MSPAEF in the presence of significant absolute mean bias would require reducing the relative contribution of the WD component. By contrast, SPAEF and E_sp already perform similarly to MSPAEF when bias values are small, but as bias increases, their lack of bias-sensitive components will limit their ability to achieve similar performance to MSPAEF.

While in this demonstration we evaluate and rank models separately for each variable and each metric, in many real-world applications, the overall model performance can be assessed using multiple variables. There are many multi-criterion model ranking techniques that can do this, such as Compromise Programming (CP) (Refaey et al., 2019; Baghel et al., 2022), Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) (Raju and Kumar, 2015) and Cooperative Game Theory (Gershon and Duckstein, 1983).

4 Conclusions

Model evaluation plays a vital role in climate model development and is crucial for understanding the limitations of these models, especially before conducting impact assessments. In this study, we present an inter-comparison exercise of several performance metrics that focus on the spatial representation of model output.

The bias-insensitive metrics, SPAEF and E_sp perform well in the cases of relatively low biases. With increasing bias, their performance generally deteriorates, something that is evident by the shift of the curve's minimum to values of standard deviation ratios greater than 1. Contrariwise, the WSPAEF metric generally has a better all-round performance than the two bias-insensitive metrics. Nonetheless, it performs poorly for skewed distributed variables as the bias increases.

To address these limitations of the existing metrics, we introduced a new metric, the Modified Spatial Efficiency (MSPAEF). MSPAEF has shown robust performance across a range of correlation, bias, and standard deviation ratio values, for variables that follow either a normal or a skewed distribution. Nevertheless, it has a slightly reduced sensitivity to changes in the standard deviation ratio near 1 when a large bias is present, particularly for normally distributed variables. This reduced sensitivity arises from the dominant influence of the large bias, primarily due to the NRMSE and relative bias components, in the final metric value.

The synthetic data examples of Sect. 3.2, demonstrate that bias-insensitive metrics are suboptimal for model evaluation when bias is a critical factor, with this limitation being particularly evident in Example 2 (Fig. 8). Conversely, for variables with large absolute values, even a small bias relative to the mean observational value can have a substantial impact on the WSPAEF metric. In contrast, the MSPAEF and the two bias-insensitive metrics are less affected, as demonstrated in Example 1 (Fig. 7). Thus, in the two examples presented in Sect. 3.2, MSPAEF was the only metric whose values consistently aligned with the intuitive assessment of model performance for both cases. This is attributed to MSPAEF being a bias-sensitive metric, unlike SPAEF and E_sp, while also accounting for the relative rather than absolute value of the mean bias, unlike the WSPAEF metric. Additionally, a distinguishing characteristic of MSPAEF is its scale-independence, meaning it is unaffected by the units of the variable, unlike the other three metrics.

In the real-world application of Sect. 3.3 using global CMIP6 models, it becomes obvious that for variables that can have large absolute value of bias, such as the total annual precipitation, the WSPAEF metric is greatly affected by it and it differs the most in the ranking compared to the other three metrics. The differences of WSPAEF with the other metrics are substantially reduced with the temperature variable, due to the smaller absolute values of the mean bias, which allows its the other components to have a larger contribution to the final value of the metric. The MSPAEF rankings are the most similar to those of E_sp for both variables, highlighting the ability of MSPAEF to evaluate spatial patterns. However, some differences in the rankings arise due to MSPAEF being a bias-sensitive metric.

Through the use of appropriate weights, the existing metrics can be adjusted to better align with the performance of MSPAEF. Nevertheless, on many occasions, such as when there is insignificant bias, their original unweighted forms often perform sufficiently well, reducing the need for such modifications.

The MSPAEF metric is not proposed as a substitute for multi-metric evaluation, which remains essential in climate model assessment. It was designed as a comprehensive similarity measure that incorporates both spatial pattern agreement and bias. It is defined in a way that emphasizes spatial pattern similarity when the relative bias is small, but it increasingly emphasizes the bias as it becomes more significant. This design allows MSPAEF to act as a balanced indicator, adapting its emphasis depending on the characteristics of the data. While this metric was developed for evaluating gridded climate model output, its underlying rationale and its inherent flexibility make it suitable for assessing the performance of other types of geoscientific or environmental models where the spatial distribution of simulated variables is expected to follow certain patterns in space. Its use alongside established measures can help provide a clearer and more consistent picture of model performance.

Appendix A: Modification of Metrics

The SPAEF, E_sp, and MSPAEF metrics were modified to take values equal to or greater than zero, for the inter-comparison of the metrics.

Specifically, the SPAEF metric was modified as follows:

$$\begin{array}{}\text{(A1)}& \text{SPAEF}=\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\gamma }-\mathrm{1}{)}^{\mathrm{2}}}\end{array}$$

The E_sp metric was modified as follows:

$$\begin{array}{}\text{(A2)}& E_{\text{sp}}=\sqrt{({\mathbf{r}}_{\mathrm{s}}-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\gamma }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}}\end{array}$$

The MSPAF metric was modified as follows:

$$\begin{array}{}\text{(A3)}& \text{MSPAEF}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{4}}}}\sqrt{(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }{)}^{\mathrm{2}}+(\mathit{\gamma }{)}^{\mathrm{2}}+(\mathit{\delta }{)}^{\mathrm{2}}}\end{array}$$

Appendix B: Probability Distributions of Precipitation from EC-Earth3 CMIP6 Historical Run

Figure B1 Probability distributions from the r1i1p1f1 variant of the historical run of the EC-Earth3 CMIP6 model. (a) Global daily precipitation for 1 January 2014. (b) Global annual precipitation for 2014. (c) Daily precipitation for the Tropics and Extra-tropics on 1 January 2014. (d) Annual precipitation for the Tropics and Extra-tropics in 2014.

Appendix C: CMIP6 models

Dix et al. (2023)Ziehn et al. (2023)Semmler et al. (2023)Xin et al. (2023)Danabasoglu (2023b)Danabasoglu (2023a)Huang (2023)Lovato and Peano (2023)Voldoire (2023b)Voldoire (2023a)Seferian (2023)Swart et al. (2023a)Swart et al. (2023b)EC-Earth Consortium (2023b)EC-Earth Consortium (2023a)YU (2023)Li (2023)Krasting et al. (2023)NASA Goddard Institute for Space Studies (2023)Ridley et al. (2023)Volodin et al. (2023a)Volodin et al. (2023b)Boucher et al. (2023)Byun et al. (2023)Hajima et al. (2023)Tatebe and Watanabe (2023)Jungclaus et al. (2023)Wieners et al. (2023)Yukimoto et al. (2023)Cao and Wang (2023)Seland et al. (2023)Bentsen et al. (2023)Tang et al. (2023)

Table C1 CMIP6 models.