When are Models Useful? Revisiting Quantification of Reality Checking

All models are wrong but some are useful (George E. P. Box ) (Box, G. E. P. Robustness in the Strategy of Scientific Model Building. Robustness in Statistics 1979, 201–236. doi:10.1016/b978-0-12-438150-6.50018-2.).

1. Introduction

The aphorism in the epigram is very popular and expresses the fact that models are only approximations of reality. The first part of aphorism, “all models are wrong” is wrong per se in a rigorous epistemological context. The meaning it purports to express could perhaps be better formulated as “models differ from reality”. They differ not only in quantitative terms—e.g., when a real value is 6 units and a model predicts 5 units. Even if a model gives a value identical to the real one (e.g. 6 units), again it differs from reality conceptually as the model is a representation (usually mathematical, simplifying and approximate) of reality—not the physical reality per se [[1]]. The model differs from the system it represents even in the case that the latter is a hardware and/or software implementation of an algorithm [[2]]. In an era where confusion has prevailed over rigour, it is essential to clarify the conceptual difference, which makes a model not “wrong” or “right”—just conceptually different from reality.

The second part of the aphorism, “some [models] are useful” is not problematic and is the subject of this paper. The usefulness of a model is an issue that needs modelling per se. In other words, it needs to be quantified by metrics that describe how good a quantified approximation of reality is. This quantification is typically based on simulation results by the model and comparison with observations. The comparison typically uses statistical or stochastic concepts such as variances and correlation coefficients. This implies that we have to take an additional modelling step, i.e. to assume that both model simulation outputs, $s$, and actual observations, $x$, are further represented as stochastic variables (or even stochastic processes), $\underset{\_}{s},\underset{\_}{x}$ (notice the notational convention to underline stochastic variables). This step may not be necessary if the model is deterministic, yet taking it allows using the advanced language and tools of stochastics, which facilitates modelling.

Standard statistical metrics of this type are the correlation coefficient, $r_{\underset{\_}{s}\underset{\_}{x}}$, between $\underset{\_}{s}$ and $\underset{\_}{x}$, and its square, $r_{\underset{\_}{s}\underset{\_}{x}}^2$, known as the coefficient of determination. These however do not provide a holistic picture of the similarity between $\underset{\_}{s}$ and $\underset{\_}{x}$ as they do not reflect the similarity (or otherwise) of the marginal distribution. For example, $\underset{\_}{s}$ and $\underset{\_}{x}$ may have a large $r_{\underset{\_}{s}\underset{\_}{x}}$, suggesting good model performance, and simultaneously a large difference in their means, suggesting poor performance because of bias.

The most common metric of a holistic type has been the Nash–Sutcliffe efficiency (NSE). It was proposed by two famous hydrologists, Nash and Sutcliffe (1970) [[3]] in a study that for a long time has been the most cited hydrological paper [[4]] (currently about 28 000 citations in Google Scholar and more than 18 000 in Scopus). Its use has been common beyond hydrology, such as in geophysics, earth sciences, atmospheric sciences, environmental sciences, statistics, engineering, data science and computational intelligence [e.g. [5]-,[6],[7],[8],[9],[10],[11],[12],[13],[14],[15],[16]]. As noted by O’Connell et al. [[17]], the discipline of hydrology has commonly been an “importer” of ideas, techniques and theories developed in other scientific disciplines. A rare exception is Hurst’s work [[18]], which was “exported” to many areas of science and technology. The model performance metric proposed by Nash and Sutcliffe (NSE) is another rare exception of “exportation”.

A different metric that has recently attracted wide attention in hydrology and beyond is the Kling–Gupta efficiency (KGE) proposed by the hydrologists Kling and Gupta [[19],[20]]. Both metrics, NSE and KGE, are expressed in terms of the first- and second-order classical moments of the variables $\underset{\_}{s},\underset{\_}{x}$ or their difference, i.e., the error $\underset{\_}{e}$. Both are dimensionless and have an upper bound, the number 1, which corresponds to perfect agreement of simulated with observed values. However, they have differences. NSE has a conceptual and rigorous definition based on the expectation of the squared error. KGE is a rather arbitrary expression, heuristically combining three indices of agreement. These also appear in NSE, if it is decomposed using stochastic algebra. It is useful that the KGE metric distinguishes the three separate indicators of agreement, but it is doubtful if their combination in one metric is useful.

Both NSE and KGE provide useful information for processes that are Gaussian or close to Gaussian, but, as will be shown in the analyses that follow, fail to perform in processes with behaviour far different from Gaussian. Therefore, we need some metrics that can provide useful information for the latter case. In addition, in non-Gaussian processes, we usually need information about the behaviour at the distribution tails. This can hardly be extracted from second-order statistical properties. On the other hand, high-order moments are unknowable if the information is extracted from data [[21]-,[22],[23],[24],[25]]. Yet, the new concept of knowable moments or K-moments [23-26] can replace the classical second-order moments, on the one hand, and extend the definition of metrics for high orders, on the other hand. This is attempted in this study, after re-examining the NSE and KGE metrics, which are based on classical moments, and locating their strengths and weaknesses (Section 2). An alternative framework based on K-moments is proposed and some synthetic examples are used to illustrate the properties of both the existing and the proposed frameworks (Section 3). In addition, an application with real-world processes, namely the precipitation process, is presented, where the real system is assumed to be described by reanalysis data for precipitation, while the models are assumed to be some popular climate models (Section 4).

2. Revisiting the Existing Framework

2.1. Nash – Sutcliffe Efficiency

The error between the simulated and actual processes is defined as

$$\underset{\_}{e} :=\underset{\_}{s}-\underset{\_}{x}$$

(1)

Based on this, the Nash – Sutcliffe efficiency is defined as

$$\mathrm{N}\mathrm{S}\mathrm{E}\mathrm{}:=1-{\displaystyle \frac{\mathrm{E}\left[{\underset{\_}{e}}^2\right]}{\mathrm{v}\mathrm{a}\mathrm{r}\left[\underset{\_}{x}\right]}}$$

(2)

The error can be decomposed as

$$\mathrm{E}\left[{\underset{\_}{e}}^2\right]=\mathrm{v}\mathrm{a}\mathrm{r}\left[\underset{\_}{e}\right]+\mathrm{E}{\left[\underset{\_}{e}\right]}^2={\sigma }_{\underset{\_}{e}}^2+{\mu }_{\underset{\_}{e}}^2$$

(3)

with ${\sigma }_{\underset{\_}{e}}^2 :=\mathrm{v}\mathrm{a}\mathrm{r}\left[\underset{\_}{e}\right],{\mu }_{\underset{\_}{e}}:=\mathrm{E}\left[\underset{\_}{e}\right]$. Hence, the NSE can be decomposed as

$$\mathrm{N}\mathrm{S}\mathrm{E}=\mathrm{E}\mathrm{V}-{\mathrm{R}\mathrm{B}}^2$$

(4)

where EV and RV are the explained variance and the relative bias, respectively:

$$\mathrm{E}\mathrm{V}=1-{\displaystyle \frac{{\sigma }_{\underset{\_}{e}}^2}{{\sigma }_{\underset{\_}{x}}^2}},\mathrm{}\mathrm{}\mathrm{R}\mathrm{B} :={\displaystyle \frac{{\mu }_{\underset{\_}{e}}}{{\sigma }_{\underset{\_}{x}}}}={\displaystyle \frac{{\mu }_{\underset{\_}{s}}-{\mu }_{\underset{\_}{x}}}{{\sigma }_{\underset{\_}{x}}}}$$

(5)

Alternatively, in the decomposition, we can substitute the statistics of $\underset{\_}{s}$ for those of $\underset{\_}{e}$ and find

$$\mathrm{E}\left[{\underset{\_}{e}}^2\right]={\sigma }_{\underset{\_}{s}}^2+{\sigma }_{\underset{\_}{x}}^2-2r_{\underset{\_}{s}\underset{\_}{x}}{\sigma }_{\underset{\_}{s}}{\sigma }_{\underset{\_}{x}}+{\mu }_{\underset{\_}{e}}^2$$

(6)

Hence:

$$\mathrm{E}\mathrm{V}=2r_{\underset{\_}{s}\underset{\_}{x}}{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}^2}{{\sigma }_{\underset{\_}{x}}^2}},\mathrm{}\mathrm{}\mathrm{N}\mathrm{S}\mathrm{E}=2r_{\underset{\_}{s}\underset{\_}{x}}{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}^2}{{\sigma }_{\underset{\_}{x}}^2}}-{\left({\displaystyle \frac{{\mu }_{\underset{\_}{s}}-{\mu }_{\underset{\_}{x}}}{{\sigma }_{\underset{\_}{x}}}}\right)}^2$$

(7)

It follows directly that, when${\mu }_{\underset{\_}{s}}={\mu }_{\underset{\_}{x}}$, the metrics $\mathrm{E}\mathrm{V}=\mathrm{N}\mathrm{S}\mathrm{E}$ are maximized for:

$${\sigma }_{\underset{\_}{s}}=\left\{\begin{array}{cc}{r}_{\underset{\_}{s}\underset{\_}{x}}{\sigma }_{\underset{\_}{x}},& r_{\underset{\_}{s}\underset{\_}{x}}\ge 0\\ 0,& r_{\underset{\_}{s}\underset{\_}{x}}\le 0\end{array}\right.$$

(8)

and their maximum value is

$${\mathrm{E}\mathrm{V}}_{\mathrm{m}\mathrm{a}\mathrm{x}}={\mathrm{N}\mathrm{S}\mathrm{E}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\left\{\begin{array}{cc}{r}_{\underset{\_}{s}\underset{\_}{x}}^2,& r_{\underset{\_}{s}\underset{\_}{x}}\ge 0\\ 0,& r_{\underset{\_}{s}\underset{\_}{x}}\le 0\end{array}\right.$$

(9)

Notably, the maximum value is not achieved when ${\sigma }_{\underset{\_}{s}}={\sigma }_{\underset{\_}{x}}$. Rather, the value that corresponds to the latter case is

$$\mathrm{E}\mathrm{V}=\mathrm{N}\mathrm{S}\mathrm{E}=2r_{\underset{\_}{s}\underset{\_}{x}}-1\le r_{\underset{\_}{s}\underset{\_}{x}}^2$$

(10)

2.1. Kling–Gupta Efficiency and Its Relationship with the Nash – Sutcliffe Efficiency

The Kling–Gupta efficiency is defined as

$$\mathrm{K}\mathrm{G}\mathrm{E}=1-\sqrt{{\left(1-r_{\underset{\_}{s}\underset{\_}{x}}\right)}^2+{\left(1-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}\right)}^2+{\left(1-{\displaystyle \frac{{\mu }_{\underset{\_}{s}}}{{\mu }_{\underset{\_}{x}}}}\right)}^2}$$

(11)

The definition is heuristic and thus KGE does not represent a formal statistic. The last term is pathological as it is often the case that ${\mu }_{\underset{\_}{x}}=0$ (e.g. when departures from the mean, known as “anomalies”, are modelled), in which KGE becomes $-\infty$. If we exclude the pathological last term, KGE becomes equivalent, but not equal, to EV in the sense that both express second-order properties.

The maximum value of KGE is achieved for ${\mu }_{\underset{\_}{s}}={\mu }_{\underset{\_}{x}},{\sigma }_{\underset{\_}{s}}={\sigma }_{\underset{\_}{x}}$, and is

$${\mathrm{K}\mathrm{G}\mathrm{E}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=r_{\underset{\_}{s}\underset{\_}{x}}$$

(12)

For ${\mu }_{\underset{\_}{s}}={\mu }_{\underset{\_}{x}},{\sigma }_{\underset{\_}{s}}\ne {\sigma }_{\underset{\_}{x}}$, after algebraic operations, we find that KGE and EV are related by

$$\mathrm{K}\mathrm{G}\mathrm{E}=1-\sqrt{{\left(1-{\displaystyle \frac{\mathrm{E}\mathrm{V}+{\sigma }_{\underset{\_}{s}}^2/{\sigma }_{\underset{\_}{x}}^2}{2{\sigma }_{\underset{\_}{s}}/{\sigma }_{\underset{\_}{x}}}}\right)}^2+{\left(1-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}\right)}^2}$$

(13)

In the limiting case that $r_{\underset{\_}{s}\underset{\_}{x}}=1$ (and ${\mu }_s={\mu }_x$), we find

$$\mathrm{K}\mathrm{G}\mathrm{E}={\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}},\mathrm{}\mathrm{}\mathrm{E}\mathrm{V}=2{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}^2}{{\sigma }_{\underset{\_}{x}}^2}}=2\mathrm{K}\mathrm{G}\mathrm{E}-{\mathrm{K}\mathrm{G}\mathrm{E}}^2\iff \mathrm{K}\mathrm{G}\mathrm{E}=1-\sqrt{1-\mathrm{E}\mathrm{V}}$$

(14)

Equations (13) and (14) are illustrated in Figure 1, where it can be seen that (a) when the explained variance EV = NSE is high, say > 0.5, the value of KGE is smaller than EV; (b) when EV = NSE is zero, the value of KGE can be as high as 0.5 (for ${\sigma }_{\underset{\_}{s}}/{\sigma }_{\underset{\_}{x}}=1)$; (c) when the EV is negative, the KGE value is less negative; (d) in general, KGE tends to decrease, in absolute value, the effectiveness metric that is provided by the EV = NSE; (e) this is further verified from the curve corresponding to ${\sigma }_{\underset{\_}{s}}/{\sigma }_{\underset{\_}{x}}$ which has a slope of 1/2, instead of 1. In brief, KGE is a less sensitive metric than NSE.

All in all, the heuristic character, the smaller sensitivity and the problematic division by the mean, which can be zero, do not favour the use of KGE and thus it will not be used further in the rest of this paper. A substitute with some mathematical meaning, namely the absolute error efficiency, AEE (rather than the squared error in NSE), is the following, computed in Appendix A:

$$\mathrm{A}\mathrm{E}\mathrm{E}\approx 1-\sqrt{{\left(1-{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}\right)}^2+2{\displaystyle \frac{{\sigma }_{\underset{\_}{s}}}{{\sigma }_{\underset{\_}{x}}}}(1-r_{\underset{\_}{s}\underset{\_}{x}})+{{\displaystyle \frac{\mathsf{\pi }}{2}}\left({\displaystyle \frac{{\mu }_{\underset{\_}{s}}-{\mu }_{\underset{\_}{x}}}{{\sigma }_{\underset{\_}{x}}}}\right)}^2}$$

(15)

Notice that the quantity $\left(1-r_{\underset{\_}{s}\underset{\_}{x}}\right)$ is not squared, as is in KGE; actually, squaring is not necessary from a mathematical point of view as it is always nonnegative. Assuming that ${\mu }_{\underset{\_}{s}}={\mu }_{\underset{\_}{x}},{\sigma }_{\underset{\_}{s}}={\sigma }_{\underset{\_}{x}}$, the resulting AEE is $1-\sqrt{2(1-r_{\underset{\_}{s}\underset{\_}{x}})}$, which takes a zero value for $r_{\underset{\_}{s}\underset{\_}{x}}=0.5$. Like in the NSE case, when${\mu }_{\underset{\_}{s}}={\mu }_{\underset{\_}{x}}$, $\mathrm{A}\mathrm{E}\mathrm{E}$ is maximized not when ${\sigma }_{\underset{\_}{s}}={\sigma }_{\underset{\_}{x}}$, but when ${\sigma }_{\underset{\_}{s}}$ satisfies equation (8).

3. Proposed Framework

3.1. A Summary of K-Moments

The methodologies discussed in Section 2 are based on first and second-order distributional properties, while the framework of classical moments cannot serve estimation of higher-order moments [22,23,26]. However, the concept of knowable moments or K-moments [23-26] can work for high -order moments.

The K-moments are defined as follows. We consider a sample of a stochastic variable $\underset{\_}{x}$, i.e., a number $p$ of independent copies of the stochastic variable $\underset{\_}{x}$, i.e., ${\underset{\_}{x}}_1,{\underset{\_}{x}}_2,\mathrm{}\dots ,\mathrm{}{\underset{\_}{x}}_p$. If we arrange the variables in ascending order, the ith smallest, denoted as ${\underset{\_}{x}}_{\left(i:p\right)}, i=1,\dots ,p$ is termed the ith order statistic. The largest (pth) order statistic is:

$${\underset{\_}{x}}_{\left(p\right)}\mathrm{}:={\underset{\_}{x}}_{\left(p:p\right)}=\max \left({\underset{\_}{x}}_1,{\underset{\_}{x}}_2,\mathrm{}\dots ,\mathrm{}{\underset{\_}{x}}_p\right)$$

(16)

and the smallest (first) is

$${\underset{\_}{x}}_{\left(1:p\right)}=\min \left({\underset{\_}{x}}_1,{\underset{\_}{x}}_2,\mathrm{}\dots ,\mathrm{}{\underset{\_}{x}}_p\right)$$

(17)

Now we define the upper knowable moment (K-moment) of order p as the expectation of the largest of the p variables ${\underset{\_}{x}}_{\left(p\right)}$:

$$K_p^{\mathrm{\text{'}}}\mathrm{}:=\mathrm{E}\left[{\underset{\_}{x}}_{\left(p\right)}\right]=\mathrm{E}\left[\max \left({\underset{\_}{x}}_1,{\underset{\_}{x}}_2,\mathrm{}\dots ,\mathrm{}{\underset{\_}{x}}_p\right)\right]$$

(18)

where $\mathrm{E}\left[\right]$ denotes expectation, and the lower knowable moment (K-moment) of order p as the expectation of the smallest of the p variables ${\underset{\_}{x}}_{\left(1:p\right)}$:

$${\overline{K}}_p^{\mathrm{\text{'}}}\mathrm{}:=\mathrm{E}\left[{\underset{\_}{x}}_{\left(1:p\right)}\right]=\mathrm{E}\left[\min \left({\underset{\_}{x}}_1,{\underset{\_}{x}}_2,\mathrm{}\dots ,\mathrm{}{\underset{\_}{x}}_p\right)\right]$$

(19)

An important property, directly resulting from their definition, is that the K-moments are ordered as follows:

$${\overline{K}}_p^{\mathrm{\text{'}}}\le \dots \le {\overline{K}}_2^{\mathrm{\text{'}}}\le {\overline{K}}_1^{\mathrm{\text{'}}}=K_1^{\mathrm{\text{'}}}=\mu \le K_2^{\mathrm{\text{'}}}\le \dots \le K_p^{\mathrm{\text{'}}}$$

(20)

These moments are noncentral and we can also define central moments as

$$K_p :=K_p^{\mathrm{\text{'}}}-K_1^{\mathrm{\text{'}}}, {\overline{K}}_p :={\overline{K}}_1^{\mathrm{\text{'}}}-{\overline{K}}_p^{\mathrm{\text{'}}}, K_p,{\overline{K}}_p\ge 0$$

(21)

As shown in [26] (chapter 6), for a stochastic variable $\underset{\_}{x}$ of continuous type, the upper K-moment of order p of $\underset{\_}{x}$, is theoretically calculated as:

$$K_p^{\mathrm{\text{'}}}=p\mathrm{E}\left[{\left(F\left(\underset{\_}{x}\right)\right)}^{p-1}\underset{\_}{x}\right]=p\underset{-\infty }{\overset{\infty }{\int }}{\left(F\left(x\right)\right)}^{p-1}x f\left(x\right)\mathrm{d}x=p\underset{0}{\overset{1}{\int }}x\left(F\right)F^{p-1}\mathrm{d}F$$

(22)

Likewise, the lower K-moment of order p is theoretically calculated as:

$${\overline{K}}_p^{\mathrm{\text{'}}}=p\mathrm{E}\left[{\left(\overline{F}\left(\underset{\_}{x}\right)\right)}^{p-1}\underset{\_}{x}\right]=p\underset{-\infty }{\overset{\infty }{\int }}{\left(\overline{F}\left(x\right)\right)}^{p-1}x f\left(x\right)\mathrm{d}x=p\underset{0}{\overset{1}{\int }}x\left(\overline{F}\right){\overline{F}}^{p-1}\mathrm{d}F$$

(23)

In these equations, $F\left(x\right)$ is the distribution function of $\underset{\_}{x}$, $\overline{F}\left(x\right) :=1-F\left(x\right)$ is its tail function and $f\left(x\right) :=\mathrm{d}F\left(x\right)/\mathrm{d}x$ is its probability density function. For discrete type variables as well as for generalizations of K-moments, the interested reader is referred to [25,26].

The unbiased estimator of the upper K-moment ${\underset{\_}{K}}_p^{\mathrm{\text{'}}}$from a sample is

$${\underset{\_}{\widehat{K}}}_p^{\mathrm{\text{'}}}=\sum _{i=1}^n{b}_{inp}\mathrm{}{\underset{\_}{x}}_{\left(i:n\right)}$$

(24)

and that of the lower K-moment is,

$${\widehat{\overline{\underset{\_}{K}}}}_p^{\mathrm{\text{'}}}=\sum _{i=1}^n{b}_{inp}\mathrm{}{\underset{\_}{x}}_{\left(n-i+1:n\right)}=\sum _{i=1}^n{b}_{n-i+1,n,p}\mathrm{}{\underset{\_}{x}}_{\left(i:n\right)}$$

(25)

where

$$b_{inp}=\left\{\begin{array}{ll}0,& i<p\\ p{\displaystyle \frac{\mathsf{\Gamma }\left(n-p+1\right)}{\mathsf{\Gamma }\left(n+1\right)}}\mathrm{}{\displaystyle \frac{\mathsf{\Gamma }\left(i\right)}{\mathsf{\Gamma }\left(i-p+1\right)}},& i\ge p\ge 0\end{array}\right.$$

(26)

where $\mathsf{\Gamma }( )$ is the gamma function. For integer moment order $p$ and $i\ge p\ge 0$ this simplifies to

$$b_{inp}=\left({\displaystyle \genfrac{}{}{0pt}{}{i-1}{p-1}}\right)/\left({\displaystyle \genfrac{}{}{0pt}{}{n}{p}}\right)$$

(27)

Based on the K-moments, we define the location (or central tendency) parameter of order $p$, $C_p$, the dispersion parameter of order $p$, $D_p$, and the dispersion ratio of order $p$, $R_p$, as:

$$C_p\mathrm{}:={\displaystyle \frac{K_p^{\mathrm{\text{'}}}+{\overline{K}}_p^{\mathrm{\text{'}}}}{2}}, D_p\mathrm{}:={\displaystyle \frac{K_p^{\mathrm{\text{'}}}-{\overline{K}}_p^{\mathrm{\text{'}}}}{2}}, R_p :={\displaystyle \frac{C_p}{D_p}}={\displaystyle \frac{K_p^{\mathrm{\text{'}}}+{\overline{K}}_p^{\mathrm{\text{'}}}}{K_p^{\mathrm{\text{'}}}-{\overline{K}}_p^{\mathrm{\text{'}}}}}$$

(28)

where $D_p\ge 0$. The least-order meaningful values thereof are:

$$C_1=K_1^{\mathrm{\text{'}}}={\overline{K}}_1^{\mathrm{\text{'}}}=\mu , D_2={\displaystyle \frac{{K_2^{\mathrm{\text{'}}}-\overline{K}}_2^{\mathrm{\text{'}}}}{2}}=K_2={\overline{K}}_2, R_2={\displaystyle \frac{K_1^{\mathrm{\text{'}}}}{K_2}}={\displaystyle \frac{{\overline{K}}_1^{\mathrm{\text{'}}}}{{\overline{K}}_2}}$$

(29)

Note that when $D_1=0,R_1=\infty .$ Also, as ${\overline{K}}_2^{\mathrm{\text{'}}}={2K}_1^{\mathrm{\text{'}}}-K_2^{\mathrm{\text{'}}}$ [26], we have $C_2=\left(1/2\right) \left(K_2^{\mathrm{\text{'}}}+{\overline{K}}_2^{\mathrm{\text{'}}}\right)=K_1^{\mathrm{\text{'}}}=C_1$.

3.2. K-moments Based Metrics of Efficiency

A perfect model will have all central tendency and dispersion parameters equal to zero: $C_p=D_p=0$ for any $p$. A good model will have nonzero values but not large. Based on the dispersion parameters, $D_p$, we can define quantities, analogous to the explained variance. We define the K-unexplained variation of order $p$ and its difference from 1, the K-explained variation of order $p$ as:

$${\mathrm{K}\mathrm{U}\mathrm{V}}_p={\displaystyle \frac{D_p\left[\underset{\_}{e}\right]}{D_p\left[\underset{\_}{x}\right]}}={\displaystyle \frac{K_p^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]-{\overline{K}}_p^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]}{K_p^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]-{\overline{K}}_p^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]}}, {\mathrm{K}\mathrm{E}\mathrm{V}}_p=1-{\mathrm{K}\mathrm{U}\mathrm{V}}_p=1-{\displaystyle \frac{D_p\left[\underset{\_}{e}\right]}{D_p\left[\underset{\_}{x}\right]}}$$

(30)

Their least-order meaningful values are

$${\mathrm{K}\mathrm{U}\mathrm{V}}_2={\displaystyle \frac{D_2\left[\underset{\_}{e}\right]}{D_2\left[\underset{\_}{x}\right]}}={\displaystyle \frac{K_2^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]-{\overline{K}}_2^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]}{K_2^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]-{\overline{K}}_2^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]}}={\displaystyle \frac{K_2\left[\underset{\_}{e}\right]}{K_2\left[\underset{\_}{x}\right]}}, {\mathrm{K}\mathrm{E}\mathrm{V}}_2=1-{\mathrm{K}\mathrm{U}\mathrm{V}}_p=1-{\displaystyle \frac{K_2\left[\underset{\_}{e}\right]}{K_2\left[\underset{\_}{x}\right]}}$$

(31)

The minimum and maximum possible values are respectively ${\mathrm{K}\mathrm{U}\mathrm{V}}_p=0$ and ${\mathrm{K}\mathrm{E}\mathrm{V}}_p=1$, and correspond to $D_p\left[\underset{\_}{e}\right]=0$. For a model that equates any $s$ with the mean of $\underset{\_}{x}$, $K_p\left[\underset{\_}{e}\right]=K_p\left[\underset{\_}{x}\right]$ for any $p\ge 2$, and ${\mathrm{K}\mathrm{U}\mathrm{V}}_p=1,{\mathrm{K}\mathrm{E}\mathrm{V}}_p=0$. Models worse than that have ${\mathrm{K}\mathrm{U}\mathrm{V}}_p$ higher than 1 and negative values of ${\mathrm{K}\mathrm{E}\mathrm{V}}_p$.

For a normal distribution, $K_2\left[\underset{\_}{e}\right]/K_2\left[\underset{\_}{x}\right]={\sigma }_{\underset{\_}{e}}/{\sigma }_{\underset{\_}{x}}$ and hence the K-moments based metrics are related to the classical explained variance by

$${\mathrm{K}\mathrm{U}\mathrm{V}}_2=\sqrt{1-\mathrm{E}\mathrm{V}}, {\mathrm{K}\mathrm{E}\mathrm{V}}_2=1-\sqrt{1-\mathrm{E}\mathrm{V}}$$

(32)

Reminding that for $r_{\underset{\_}{s}\underset{\_}{x}}=1$, $\mathrm{K}\mathrm{G}\mathrm{E}=1-\sqrt{1-\mathrm{E}\mathrm{V}}$ (see Figure 1), we notice the identity of ${\mathrm{K}\mathrm{E}\mathrm{V}}_2$ and $\mathrm{K}\mathrm{G}\mathrm{E}$ (as functions of EV) for this case and further observe that Equation (32) holds for any $r_{\underset{\_}{s}\underset{\_}{x}}$ for normal distribution, while that for KGE only holds for $r_{\underset{\_}{s}\underset{\_}{x}}=1$.

The bias is a separate characteristic and it would be better treated based on a different statistic. The ratio $R_p$ is an appropriate metric for it. Alternatively—and in an analogous manner to ${\mathrm{K}\mathrm{U}\mathrm{V}}_p$—we can define the K-bias of order p as

$${\mathrm{K}\mathrm{B}}_p={\displaystyle \frac{K_p^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]+{\overline{K}}_p^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]}{K_p^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]-{\overline{K}}_p^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]}}$$

(33)

with a special case for $p=2$

$${\mathrm{K}\mathrm{B}}_2={\displaystyle \frac{K_2^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]+{\overline{K}}_2^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]}{K_2^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]-{\overline{K}}_2^{\mathrm{\text{'}}}\left[\underset{\_}{x}\right]}}={\displaystyle \frac{K_1^{\mathrm{\text{'}}}\left[\underset{\_}{e}\right]}{K_2\left[\underset{\_}{x}\right]}}$$

(34)

It appears more natural for the two quantities describing dispersion and bias to be thought of as small as possible in order for a model to be regarded as good. To fulfil this desideratum, the metrics of choice would be the ${\mathrm{K}\mathrm{U}\mathrm{V}}_p$ and ${\mathrm{K}\mathrm{B}}_p$. However, as the aim of this paper is to provide metrics with behaviour similar to the existing ones, in what follows we use ${\mathrm{K}\mathrm{E}\mathrm{V}}_p$ rather than the more natural ${\mathrm{K}\mathrm{U}\mathrm{V}}_p$. In the case that the unexplained variation and bias are to be combined, a relevant expression, called the K-moment-based absolute error efficiency and derived in Appendix A, is:

$$\mathrm{K}\mathrm{A}\mathrm{E}\mathrm{E}\approx 1-\sqrt{{\left({\displaystyle \frac{K_2\left[\underset{\_}{e}\right]}{K_2\left[\underset{\_}{x}\right]}}\right)}^2+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{K_1^{\text{'}}\left[\underset{\_}{e}\right]}{K_2\left[\underset{\_}{x}\right]}}\right)}^2}=1-\sqrt{{\mathrm{K}\mathrm{U}\mathrm{V}}_2^2+{\displaystyle \frac{1}{2}}{\mathrm{K}\mathrm{B}}_2^2}$$

(35)

3.2. Possible Transformations of Data

Figure 2 (upper) compares two synthetic time series consisting of 1000 data values, an original $x_i$ and a simulated $s_i$. These synthetic series were constructed as follows. First a series $v_i$ was generated from the Hurst-Kolmogorov model with Hurst parameter 0.95 and normal distribution N(0,1). Then the series $v_i$ was smoothed by a linear filter with a triangular shape, with a peak value of 1 and zero values at times 10 and 20 time steps before and after the time of the peak, respectively, thus producing a series $y_i$. Subsequently, a series $z_i$ was generated by adding to $y_i$ a series $u_i$ again generated from the Hurst-Kolmogorov model with Hurst parameter 0.95 and normal distribution N(0,$\sqrt{2}{\sigma }_y$), where ${\sigma }_y$ is the standard deviation of $y_i$. Finally, the two latter series were exponentiated, $x_i=\exp (a{z}_i-{\mu }_z),s_i=\exp (a{y}_i-{\mu }_y)$, with $a=0.3$, and taken as original and simulated series after rounding to one decimal point. Both $x_i$ and $s_i$ exhibit long-range dependence and are lognormally distributed. The former is rough, due to the component $u_i$, while the latter is smooth.

From the visual depiction of Figure 2 (upper) it turns out that the model performance is very poor. This is also reflected in all metrics shown in Table 1 (first row). However, the poor performance is mostly due to the log-normal distribution of $x_i$, which yields frequent high peaks. The simulated series $s_i$ does not capture these peaks. If we take the logarithmic transformations of both $x_i$ and $s_i$, then, as seen in Figure 2 (middle), there is some resemblance of the original and simulated series, which is also reflected in the metrics of Table 1 for the log-transformed series.

This fact suggests that if the series examined are not normally distributed, a transformation of the data may be necessary before assessing the suitability of the model through appropriate indices. The logarithmic transformation may not be appropriate for all cases and also has a problem, namely the fact that it diverges to minus infinity when the original value is zero. The latter problem was resolved in an ad hoc manner to calculate the values in Table 1, i.e. by mapping the zero values in $x_i$ and $s_i$ to the logarithm of 0.01 (one order of magnitude lower than the minimum value, which is 0.1).

A proper transformation that remedies these problems, by being general and avoiding ad hoc considerations, is the following [26] (section 2.10), which we denote as lambda (λ) transformation:

$$x_i^{*} :=\lambda \ln \left(1+x_i/\lambda \right)$$

(36)

and likewise for $s_i$. For low values of $x_i\ll \lambda$, including $x_i=0$, this maps $x_i$ to itself, while for large $x_i$ it maps it to a linear function of $\ln x_i$. The parameter λ is assumed to be the same for both $x_i$ and $s_i$. We can estimate it numerically by maximizing one of the model efficiency metrics.

Furthermore, for the simulated series, unless the model is physically based and its simulation results have some physical meaning, we may give it some more degrees of freedom by the additional parameters α and β using the transformation

$$s_i^{**}:=\alpha +\beta \lambda \ln \left(1+s_i/\lambda \right)$$

with default values $\alpha =0,\beta =1$. Again, these are obtained by optimization together with the optimization of $\lambda$, noting that $\lambda$ applies to both $x_i$ and $s_i$, but α and β apply to $s_i$ only. Table 1 gives the optimized values of $\lambda$ for the default values of α and β, as well as the optimized values of all three parameters, and the resulting optimal metrics of model efficiency. Figure 2 (lower) shows the transformed time series $x_i^{*}$ and $s_i^{**}$, where a good agreement between the two is seen, in contrast to the upper panel of the same figure.

The metrics optimized for the two λ-transformed series of Table 1 are the KEV₂ when the default default values $\alpha =0,\beta =1$ are used and the KAEE otherwise. These are not the only options as high-order metrics could also have been chosen to be optimized. The higher-order metrics when the KEV₂ is optimized are shown in Figure 3 (left), also in comparison to those of the untransformed series. As seen in the figure, the performance deteriorates with the increase of the moment order. On the other hand, Figure 3 (right) shows that the performance, in terms of the KEV₂ metric, improves when the time scale increases.

Figure 3 does not give any information about the model bias. This is provided in Figure 4 in terms of both the ratio$R_p$ and the K-bias${\mathrm{K}\mathrm{B}}_p$. As seen in the graphs, both indices practically provide the same information and there is no substantial variation of the bias metrics with order p. However, the λ transformation reduces the bias of the original series substantially.

4. Real World Case Study

To present a large-scale case study of hydrological interest, we use the results of climate models for precipitation, which have been very popular and widely used in so-called climate impact studies, but without proper testing whether they are useful or not. The climate models (also known as global circulation models—GCM) used belong to the last-generation Coupled Model Intercomparison Project (CMIP6) and their outputs for precipitation were retrieved from the Koninklijk Nederlands Meteorologisch Instituut (KNMI) Climate Explorer [[26],[27]]. The outputs from the 37 models listed in Table 2 were available on a monthly scale and were aggregated to annual and over-annual scales.

As a time series representing reality, the gridded data of the ERA5 reanalysis were used [[28]]. This is the fifth-generation atmospheric reanalysis of the European Centre for Medium Range Weather Forecasts (ECMWF), where the name ERA refers to ECMWF ReAnalysis. ERA5 has been produced as an operational service and its fields compare well with the ECMWF operational analyses. It combines vast amounts of historical observations into global estimates using advanced modelling and data assimilation systems. The data are available for the period 1940-now at a spatial resolution of 0.5° globally and were retrieved using the Web-based Reanalyses Intercomparison Tools (WRIT) [[29]] made available by the USA National Oceanic and Atmospheric Administration (NOAA). The comparisons of models and reality, represented by ERA5, were made for the period 1940-2023 (84 years), separately for the North Hemisphere (NH) and the South Hemisphere (SH). A visual comparison of the time series is presented using spaghetti graphs in Figure 5 on the annual scale and Figure 6 on 8-year scale. The latter was selected as the maximum climatic scale that allows 10 data points, so that statistics can be estimated.

A prominent characteristic seen in the spaghetti graphs is the large bias of models, which is mostly negative for NH and mostly positive for SH. Different models have largely different bias. The large bias in precipitation certainly reflects inappropriate modelling of the physical processes related to the hydrological cycle, starting with latent heat and evaporation.

Nonetheless, Figure 7 shows that on a hemispheric basis there is correlation between models and reality, with an average of 0.31 for NH and 0.11 for SH. An interesting property is that each model’s precipitation at NH is negatively correlated to that of the same model for the SH, with an average correlation of –0.61. This model property, however, does not correspond to reality: if this correlation is estimated from the ERA5 data it is practically zero (–0.03).

These correlations are not enough to suggest usefulness of the models in terms of the explained variance. As seen in Figure 8, both classical and K-explained variance are mostly negative on the annual scale. Yet, positive values appear on the 8-year scale. However, if we also consider the bias, which is shown in Figure 9, the total efficiency metrics, NSE and KAEE take highly negative values, which prohibit any usefulness of the climate models for hydrological purposes.

A most relevant question is whether or not some of the models have relatively good performance in general. To study this question we first assess which of the models have the best performance, in a Pareto optimality sense, for both hemispheres. To this aim, we plot in Figure 10 the cross-performance of the GCMs at both hemispheres, in terms of the K-explained variation, KEV₂, on annual and 8-year time scales. On the annual scale, the models show mostly negative explained variation, that is, poor performance. On the 8-year scale, the performance is improved. The two models with the best performance appear to be the CMCC-CM2-SR5 and FGOALS-f3-L.

The change of the performance with moment order and time scale for model CMCC-CM2-SR5 is shown in Figure 11, where it is seen that for small time scales the performance is not good in none of the hemispheres. Figure 12 shows that the performance at the annual scale can be slightly improved by applying the transformation of Equation (37), but this is accompanied by a worsening of the performance at large time scales. Interestingly, the improvement at the annual scale is due to the linear part of the transformation, namely on the parameter $\beta =0.41\ne 1$, and not due to the logarithmic part.

The next step is to choose an area smaller than an entire hemisphere and assess the performance of the two “best performing” models in this area. Given that ERA5 is developed in Europe and hence expected to be more accurate in this area, we chose a spherical rectangle that contains Europe, namely that defined by the coordinates 11° W, 40° E, 34° N, 71° N. The time series of the two models in question, integrated over this area, are shown in Figure 13, in comparison with the ERA5 time series. The visual comparison is not encouraging in terms of the agreement of models with reality.

Even without considering the bias, which is substantial, i.e., by only considering the explained variation, the results are rather disappointing, with ${\mathrm{K}\mathrm{E}\mathrm{V}}_2^{\left(\kappa \right)}$ not exceeding 0.1 at any scale $\kappa$.

5. Discussion and Conclusions

The classical Nash–Sutcliffe efficiency appears to be a good metric of the appropriateness of a model. Yet its fusion of two different characteristics, the explained variance and the bias, is not always an advantage. The bias could be a very important characteristic if a model is physically based and the bias reflects a violation of a physical law (e.g. conservation of mass of energy). In such cases, a large bias would suffice to reject a model, even if it captures the variation patterns.

In other cases, in which the model is of conceptual mathematical, rather than physical, type, the bias can be easily removed by a shift of the origin. In such cases, a nonlinear transformation of the observed and modelled series, accompanied by a linear transformation of the simulated series (Equations (36) – (37)), can potentially improve the agreement between the model and reality. It is suggested that in such cases, the quantified assessment of model usefulness be based on the metrics of the transformed series.

The typical metrics that are currently used to assess model performance are based on classical statistics up to second order. This is not a problem when the processes are Gaussian, but most hydrological processes are non-Gaussian. The concept of knowable moments (K-moments) offers the basis to extend the performance metrics to high orders, up to the sample size. The two metrics proposed, the K-unexplained variation, ${\mathrm{K}\mathrm{U}\mathrm{V}}_p$, and the K-bias, ${\mathrm{K}\mathrm{B}}_p$, both based on K-moments of the model error, provide ideal means to assess the agreement of models with reality; the closer to zero they are, the better the agreement. The lowest order on which they are evaluated is $p=2$, which represents second-order properties, but also using higher orders gives useful information on the agreement of the entire distribution functions.

The real-world application presented is a large-scale comparison of climatic model outputs for precipitation with reality for the last 84 years. It turns out that the precipitation simulated by the climate model does not agree with reality on the annual scale but there is some improvement on larger time scales at a hemispheric basis. However, when the areal scale is decreased from hemispheric to continental, i.e. when Europe is examined, the model performance is poor even at large time scales. Therefore the usefulness of climate model results for hydrological purposes is doubtful.