3.1 A cautionary note on resampling

One difficulty when determining the required ensemble size for a specific question is the chosen sampling approach: in this study, we generate synthetic ensembles of different ensemble sizes by randomly sampling members from a 200-member ensemble without replacement. Samples generated in this way are not fully independent when approaching the full ensemble size. For example, two random samples of 190 out of the available 200 members will share most of their members. This resampling introduces a problem when the signal is defined by using the full ensemble. Any subsample that is close to the full ensemble size will then indicate that the ensemble size is sufficient by construction.

The resampling problem occurs with any limited sample. At some point, the 1000 random subsamples are not independent anymore because they share many of the randomly drawn members from the full ensemble. Therefore, they look more similar not only to each other but also to the 200-member mean. To demonstrate how this resampling affects our estimate of the error, we deliberately reduce the size of the ensemble, for instance by only using the first 150 members and repeating the analysis. In an empirical analysis, we find that samples using more than 50 % of the available ensemble size to generate random samples lead to a substantial bias in the error estimate. We therefore recommend treating results indicating that, for example, more than 100 out of 200 members are required with caution because the true required ensemble size might be much larger. A more detailed discussion is provided in Appendix A.

3.2 A recipe for estimating ensemble size

Based on the example introduced in this section, we suggest the following approach to derive a robust estimate of the required ensemble size for any application. This method can be applied either to one of the existing large ensembles, as shown above for the MPI-GE, or to a long control run, which is available for all models participating in CMIP. We summarise the method in five steps before applying it to several examples in the next section:

1. Define the question to be addressed (isolate a forced response, quantify variability, or detect a change in variability).

2. Choose an error metric (e.g. RMSE or variance across samples) and an upper threshold based on the maximum error that is acceptable in the specific application.

3. Estimate the error for different ensemble sizes by subsampling a long control run or a large ensemble of transient simulations.

4. Determine the minimum ensemble size that is required to reduce the error below the threshold chosen in step 2.

5. If the ensemble size determined in this way is less than 50 % of the available sample size (e.g. 50 members when subsampling a 100-member ensemble), then the estimated required ensemble size provides a robust estimate for the specific question and model investigated. If the estimated required ensemble size is larger than 50 % of the available sample size, then the estimate is biased low and the true required ensemble size could be substantially larger.

4.1 Quantifying the forced response

The forced response shown in Fig. 1 contains various signals. The most prominent signal is the long-term warming trend caused by anthropogenic greenhouse gas emissions. On shorter timescales, volcanic eruptions lead to a cooling of the global mean surface temperature.

In the first example, we continue to use the RMSE to quantify how well the entire forced response is estimated, but we move from the global mean to the regional forced response in near-surface air temperature in the historical runs from the MPI-GE. In Fig. 3a–e, the expected RMSE for each grid point is shown for ensemble sizes of 3, 5, 10, 50, and 100 members. This analysis is equivalent to the computation of the mean RMSE for GSAT (black line in Fig. 2), but applied to each grid point separately. The RMSE is computed as the mean difference between 100 samples and the 200-member mean. When the ensemble mean is based on just three members, the expected error in the estimated forced response is large over land regions, in particular in the Northern Hemisphere. Over the ocean, the RMSE is already small in many regions. Increasing the ensemble size reduces the error. At 50 members, the error is small in most regions of the globe. Because 50 members is smaller than 50 % of the maximum ensemble size (200 members), the error estimate for this ensemble size is reliable.

Figure 3(a–e) The mean RMSE for the forced response in historical monthly mean near-surface air temperature of MPI-GE for (a) 3, (b) 5, (c) 10, (d) 50, and (e) 100 ensemble members relative to the 200-member mean, globally. The RMSE shown here is the mean from 100 random samples without replacement. (f–j) Required ensemble size to capture the 200-member mean forced response in historical monthly mean near-surface air temperature dependent on the acceptable error of (f) 0.1, (g) 0.2, (h) 0.3, (i) 0.5, and (j) 1.0 ^∘C.

To estimate how many members are sufficient to reduce the error below a critical threshold, we first need to determine what is an acceptable error as outlined in step 2 of the recipe. This choice will depend on the region of interest and the accuracy with which the forced response needs to be quantified. In Fig. 3f–j, we show how many members are necessary to estimate the forced response in near-surface air temperature for five acceptable errors that were chosen for illustrative purposes. If the acceptable error (RMSE) is 0.1 ^∘C, 10–30 ensemble members are sufficient over the tropical ocean, while more than 50 ensemble members are required over most land regions. Beyond 100 members, the resampling problem inhibits reliable estimates of the sufficient ensemble size. For an acceptable error of 0.25 ^∘C, fewer than 10 members are sufficient over most ocean regions, while more than 50 members are required over high-northern-latitude land regions. For an acceptable error of 0.5 ^∘C, only high-latitude land regions require a large ensemble, while the forced response over ocean and land regions at lower latitudes can be estimated with fewer than 10 members.

Conversely for rainfall, the error in estimating the forced signal when using a small ensemble is larger over the tropics than over the higher latitudes (Fig. 4a–e). The largest errors can be found over the Indian Ocean and western tropical Pacific. Similar to temperature, a 50-member ensemble shows very small errors across the globe.

In Fig. 4f–j we show how many members are necessary to estimate the forced response with an acceptable error of 0.1, 0.2, 0.3, 0.5, and 1 mm d⁻¹. For an acceptable error of 0.2 mm d⁻¹, some ocean regions require more than 100 members to capture the forced rainfall response with the required accuracy, while fewer than 20 members are sufficient over northern Africa and Eurasia. Over large parts of North America, between 20 and 40 members are required to estimate the forced rainfall response. For an acceptable error of 0.5 mm d⁻¹, 20 to 40 members are required over the Indian Ocean and western tropical Pacific, while fewer than 10 members are sufficient elsewhere.

Figure 4(a–e) The mean RMSE for the forced response in historical monthly mean total precipitation of MPI-GE for (a) 3, (b) 5, (c) 10, (d) 50, and (e) 100 ensemble members relative to the 200-member mean, globally. The RMSE shown here is the mean from 100 random samples without replacement. (f–j) Required ensemble size to capture the 200-member mean forced response in historical monthly mean total precipitation dependent on the acceptable error of (f) 0.1, (g) 0.2, (h) 0.3, (i) 0.5, and (j) 1.0 mm d⁻¹.

For the example in Figs. 3 and 4, the objective was to isolate the full forced response in a time series, defined as the 200-member ensemble mean time series at every grid point. The full forced response includes all external forcings, both natural and anthropogenic. In many applications, the objective might be to isolate a specific feature of the forced response rather than all components. In the following two examples, we will demonstrate how to estimate the required ensemble size needed to isolate the global warming trend in the 20th century and the global cooling after a major volcanic eruption.

The global warming signal follows a much simpler trajectory than the forced response to all external forcings (cf. Fig. 1). Here, we fit a linear trend to the historical time series for 1920 to 2005 and define the 200-member mean as the true forced warming trend. Over the 68-year period from 1920 to 2005, the model warms by 0.65 K (Fig. 5). We acknowledge that a linear trend may not represent the anthropogenic warming accurately but use this definition to illustrate how a specific aspect of the forced response can be investigated.

Figure 5Linear warming trend from 1920 to 2005 for different ensemble sizes shown as a linear trend fitted to the ensemble mean. Black lines show maximum and minimum 86-year ensemble mean temperature trend from 1000 random samples. Errors are shown as percentage of the 200-member ensemble mean temperature trend.

Download

We subsample the ensemble for smaller ensemble sizes to generate forced warming trends for smaller ensemble sizes. While the trends in a single realisation can be anywhere in the range from 0.4 K to more than 0.8 K warming over 68 years, increasing the ensemble size to even five members leads to a significant reduction in the error (Fig. 5). The warming trend in every 10-member ensemble is within the 20 % range (±10 %, cyan dashed lines) of the true warming trend, indicating that ensembles with 5–10 members can provide a good estimate of the forced linear warming trend. While an error within the 20 % range of the true signal may be sufficient for some applications, the acceptable error for other applications might be larger or smaller and result in a smaller or larger acceptable ensemble size. For an acceptable error of ±15 %, five ensemble members would be sufficient, while for an acceptable error of ±5 % at least 25 ensemble members are required. All of these error estimates are below 100 members and therefore not dominated by the resampling problem.

Figure 6GSAT cooling after Krakatoa eruption for different ensemble sizes shown as the ensemble mean temperature difference between 1882 and 1884. Black lines show maximum and minimum temperature response from 1000 random samples. Errors are shown as percentage of the 200-member ensemble mean temperature response.

Download

For signals on shorter timescales, the required ensemble size can be quite different. In Fig. 6 we analyse the GSAT cooling after the Krakatoa eruption in 1883. The forced cooling is quantified as the difference between 1884, the year after the eruption, and 1882, the year before the eruption. The 200-member mean shows a forced cooling of −0.34 K after the eruption. Due to internal variability, a single realisation can even show a warming after the volcanic eruption. More than one member is required for the ensemble mean to capture a cooling in all samples. However, the ensemble mean cooling for five members can still exceed the range from −0.2 to −0.5 K. More than 50 ensemble members are necessary to estimate the forced cooling within ±15 % of the true forced cooling, and approximately 100 members are required to reduce the error below ±10 %. Due to the resampling problem, we cannot derive a robust estimate for the ensemble size required to reduce the error to less than ±5 %. While the analysis in Fig. 6 suggests that 150 members would be sufficient for a ±5 % error, this number is close to the full ensemble size of 200 members and therefore biased low. The true required ensemble size to reduce the error to ±5 % is likely larger than 150 members.

These examples demonstrate that the required sample size to estimate the forced response depends on the region and variable (Figs. 3 and 4), as well as the feature of interest in the forced response (Figs. 5 and 6). Whereas for some applications five members are sufficient to reduce the error to an acceptable magnitude, other applications require at least 50 members. A robust estimate for the forced response is given by the ensemble mean when averaging over the ensemble attenuates internal variability sufficiently (Frankcombe et al., 2018). The number of members required for this depends both on the magnitude of the forced signal and the magnitude of internal variability, as well as on the acceptable error for a specific application.

4.2 Quantifying internal variability

While quantifying the forced response only requires a robust estimate of the mean, quantifying internal variability requires more members because higher-order moments of the distribution need to be estimated. In the following two examples, we use the second statistical moment of the distribution, the standard deviation, to quantify internal variability. We note that, if the distribution deviates from a normal distribution, only using the standard deviation to quantify internal variability may not be sufficient.

Here, we investigate internal variability in two regions: the tropical Pacific, where the variability is primarily driven by ENSO, and the central United States (34–46^∘ N, 116–96^∘ W). The tropical Pacific region shows substantial variability on interannual to decadal timescales. Previous work has demonstrated that large sample sizes are necessary to quantify ENSO variability (Maher et al., 2018; Wittenberg, 2009). As a second region, we analyse temperature variability over the central United States. We hypothesise that these two regions should have different requirements for the ensemble size, with a smaller required ensemble size for the central United States than the tropical Pacific to stay within an acceptable error range.

For the following examples we use the 2000-year pre-industrial control simulation from the MPI-GE. The advantage of this approach, in contrast to the examples for the forced response, is that the required ensemble size can be estimated for any model without needing a large ensemble to be available. The disadvantage is that, when using the pre-industrial control simulation, we assume that internal variability does not change under global warming.

We quantify ENSO variability by using the December, January, and February (DJF) variability in the Niño3.4 box (5^∘ N–5^∘ S, 170–120^∘ W). To ensure that ENSO variability on interannual to multi-decadal timescales is sampled, we use the Niño3.4 standard deviation for a 100-year period. The standard deviation, as computed for the full 2000-year time series, is used as the truth in this context and indicated by the horizontal black line in Fig. 7a. To generate synthetic ensemble members, we split the pre-industrial control simulation into overlapping 100-year segments. Each segment is used as one ensemble member, and the temporal standard deviation over the 100-year segment represents ENSO variability for this member. For an ensemble size of one, the spread in ENSO variability seen in Fig. 7a indicates that individual 100-year periods can have substantially more or less variability than the reference value based on the full control run.

Figure 7We show for increasing ensemble sizes the (a) ENSO variability in the Niño3.4 box (5^∘ N–5^∘ S, 170–120^∘ W) calculated over 100-year periods, (b) central United States variability (34–46^∘ N, 116–96^∘ W) calculated over 100-year periods, (c) ENSO variability in the Niño3.4 box calculated over 30-year periods, and (d) central United States variability calculated over 30-year periods. All indices are calculated from the 2000-year MPI-GE control run. Each index is calculated as a running value at each time step in the control. ENSO indices are calculated for DJF, and United States indices are calculated for the annual mean. Ensembles of 1 to 120 members are created by randomly sampling the control simulation without replacement. For each ensemble size we create 1000 artificial ensembles. The estimated true value is calculated by using the entire 2000 years of the control and is shown as a horizontal black line. The maximum and minimum values of each index from the 1000 samples are shown as solid black lines. Varying error thresholds are shown as horizontal coloured lines.

Download

To account for this centennial modulation of ENSO variability, the ENSO variability in multiple ensemble members can be averaged to get a more accurate estimate of the average ENSO variability. We simulate different ensemble sizes by averaging over randomly chosen members for a given ensemble size and repeat this 1000 times. By using a five-member mean, the error of the estimated variability in all samples is within ±15 % of the true value. To reduce the error below ±10 %, 10 ensemble members are sufficient. To improve the accuracy so that the ENSO variability estimate is within ±5 % of the truth, nearly 50 ensemble members are necessary.

For a region with less variability, much smaller ensemble sizes are sufficient to obtain a similar accuracy. For annual mean central US temperatures (Fig. 7b) any individual realisation is within ±15 % of the truth and 10 members are sufficient to increase the accuracy to the ±5 % range around the truth, whereas 50 members are necessary for ENSO. This emphasises that for some regions and quantities a moderate ensemble size or even a single realisation can be sufficient to quantify internal variability.

In both examples, the long sampling period of 100 years increases the sample size and thereby improves the accuracy for individual realisations. This is useful if the objective is to quantify variability when stationarity can be assumed, but it can be problematic if the objective is to identify a change in variability, such as changes in ENSO characteristics under global warming. A more detailed discussion of estimating ENSO variability, in particular using the ensemble dimension instead of the time dimension in transient simulations to quantify internal variability, can be found in Maher et al. (2018) and Haszpra et al. (2020).

4.3 Quantifying changes in internal variability

To quantify changes in internal variability, we need a robust estimate of internal variability both for a reference period and for a period where we want to investigate a potential change in variability (e.g. a pre-industrial control state and a time period in a future scenario). This problem is more challenging than the previous examples because the errors for the variability estimates of the two time periods add up. To demonstrate this, we use the internal variability of September Arctic sea ice area as an example. Previous work has shown that the internal variability in Arctic sea ice area first increases under warming, before it approaches zero when most of the Arctic sea ice has melted (Goosse et al., 2009; Olonscheck and Notz, 2017). We analyse the 100 members from the 1 % CO₂ scenario from the MPI-GE and use the ensemble standard deviation as an estimator of internal variability. After 120 years, nearly all ensemble members show a completely ice-free Arctic in September (Fig. B1a). The internal variability increases from model year 1 to year 80, before it sharply drops, reaching zero around year 120 (Fig. B1b).

Figure 8Change in internal variability of September Arctic sea ice is from the first decade to years 71–80 in a 1 % CO₂ experiment. For different ensemble sizes, we compute the ensemble standard deviation and then average for the first decade and years 71–80 before computing the difference. Black lines show maximum and minimum change in variability from 1000 random samples. Errors are shown as percentage of the 100-member variability change.

Download

Here we focus on the increase in variability from the beginning of the simulation to year 80 and ask how many ensemble members are necessary to robustly quantify this change in internal variability. To increase the sample size, we use a decadal mean of the ensemble standard deviation rather than a single year. We then compute the difference in internal variability between the two time periods for ensemble sizes between 3 and 100 members. Figure 8 shows the range of this change in internal variability from 1000 random samples. To quantify the change in variability within ±15 % of the true value (here defined as the internal variability change estimated with 100 members), 50 ensemble members are necessary. An error of less than ±10 % and ±5 % is only reached beyond 50 members. Due to the effect of resampling beyond 50 members, we cannot estimate the required ensemble size for these error thresholds from the 100-member ensemble used here. For very small ensemble sizes, the estimate of the variability change may even show the opposite sign of the true change, i.e. a decrease in internal variability.

The large number of ensemble members required to robustly quantify this change in variability shows that identifying a change in internal variability requires the largest ensemble size of all examples shown in this study, even when using decadal averaging to increase the sample size. This is because a robust estimate of a change in internal variability requires a clean separation of internal variability from the forced response and a robust estimate of internal variability for two different time periods. Errors in any of these estimates will propagate to the estimated change in variability, thereby making it more challenging. A small forced change in internal variability will further complicate this analysis.

A first estimate for the magnitude of a detectable change in internal variability can be derived from the control run (as in Fig. 7). Any change in variability that is smaller than the uncertainty of the estimated internal variability for a given ensemble size is not detectable. We note that this method can also be used to add error bars to estimates of forced changes in internal variability under climate change in small ensembles or single realisations from CMIP and hence determine the robustness of results.