Evidence for climate change in the satellite cloud record

Clouds substantially affect Earth’s energy budget by reflecting solar radiation back to space and by restricting emission of thermal radiation to space. They are perhaps the largest uncertainty in our understanding of climate change, owing to disagreement among climate models and observational datasets over what cloud changes have occurred during recent decades and will occur in response to global warming. This is because observational systems originally designed for monitoring weather have lacked sufficient stability to detect cloud changes reliably over decades unless they have been corrected to remove artefacts. Here we show that several independent, empirically corrected satellite records exhibit large-scale patterns of cloud change between the 1980s and the 2000s that are similar to those produced by model simulations of climate with recent historical external radiative forcing. Observed and simulated cloud change patterns are consistent with poleward retreat of mid-latitude storm tracks, expansion of subtropical dry zones, and increasing height of the highest cloud tops at all latitudes. The primary drivers of these cloud changes appear to be increasing greenhouse gas concentrations and a recovery from volcanic radiative cooling. These results indicate that the cloud changes most consistently predicted by global climate models are currently occurring in nature.

2002-2014 CERES albedo and the 1985-1989 ERBS albedo. All observational records agree that cloud amount and albedo increased over the northwest Indian Ocean, the northwest and southwest tropical Pacific Ocean, and north of the Equator in the Pacific and Atlantic oceans. Cloud amount and albedo decreased over mid-latitude oceans in both hemispheres (especially over the North Atlantic), over the southeast Indian Ocean, and in a northwest-to-southeast line stretching across the central tropical South Pacific. MAC-LWP exhibits a similar trend pattern in liquid water path during the period 1988-2014 (Extended Data Fig. 1a). Shifting the start or end time of trend calculation by several years has little impact on the spatial pattern of change. Similar patterns occur for differences in total cloud amount and albedo between the periods 1985-1989 and 2003-2009, during which ISCCP, PATMOS-x, and ERBS/CERES completely overlap (Extended Data Fig. 2).
Are the observed cloud changes solely a manifestation of natural internal variability or are they also a response to external radiative forcing of the climate system? We addressed this question by examining simulations from the Coupled Model Intercomparison Project Phase 5 (CMIP5) multi-model dataset 17 . Historical simulations included anthropogenic greenhouse gas concentrations, ozone, land-use changes, anthropogenic aerosols, volcanic aerosols, and solar output and thus represent our best estimate of the climate response to recent external radiative forcing (Extended Data Table 1). Figure 1c displays the spatial distribution of trends in ensemble mean modelled total cloud amount during the 27-year period 1983-2009 for all radiative forcings (ALL). Observations and models exhibit widespread agreement on which areas have increasing and decreasing cloud amount (Fig. 1d). Table 1 lists the spatial correlation between observed and simulated cloud trends.
Could natural internal variability alone produce the correlation between the observed and simulated cloud trend patterns? We assessed the likelihood of this outcome by examining cloud trends during 27-year periods from CMIP5 preindustrial simulations without external radiative forcing (Extended Data Table 2). Calculating the spatial correlation between the ensemble mean ALL trend pattern and the trend pattern of each 27-year preindustrial period generates a probability distribution of correlation values arising purely from natural internal variability (Extended Data Fig. 3). We found that no 27-year period in more than 15,000 years of preindustrial simulations exhibits a correlation coefficient as positive as that between the observed and ensemble mean ALL trend patterns, suggesting that external radiative forcing was a driving factor in large-scale cloud changes from the 1980s to the 2000s.
One prominent feature of Fig. 1, and a robust prediction by climate models, is the widespread reduction in cloudiness at midlatitudes 2,10,12,18 . Figure 2a, b shows trends in zonal mean total cloud amount during the period 1983-2009 for ISCCP and PATMOS-x, and Fig. 2c shows zonal mean differences between the 2002-2014 CERES albedo and the 1985-1989 ERBS albedo. Every observational record exhibits a decline in cloud amount or albedo at mid-latitudes in both hemispheres that is nearly always statistically significant. The oceanonly MAC-LWP dataset also reports less liquid water path around 40° N and 40° S (Extended Data Fig. 1b). Previous research found evidence for tropical expansion in recent decades 19 . Reduced cloudiness around 40° N and 40° S is consistent with a poleward expansion of the subtropical dry zone cloud minimum and poleward retreat of the storm-track cloud maximum. Figure 2d displays trends in zonal mean total cloud amount during the period 1983-2009 from the ALL simulations. Most individual simulations exhibit reduced cloud amount in the mid-latitudes of both hemispheres, and the ensemble mean trends are statistically significant (P < 0.05 two-sided). Furthermore, the majority of simulations reproduce the observed increase in cloud amount and albedo occurring in the northern tropics. The spatial correlation between observed and simulated zonal cloud trends is highly significant (Table 1 and Extended Data Fig. 3).
Since the correction procedures applied to the satellite datasets removed any real global mean change that might be present, for maximum comparability we subtracted the 60° S-60° N average change in total cloud amount from the model output before creating Figs 1c and d, and 2d. Without this adjustment, the ALL ensemble mean cloud amount averaged over 60° S-60° N decreases by 0.13% over 25 years. Although highly statistically significant (P < 0.0001 two-sided), the modelled reduction in 60° S-60° N average cloud amount during the period 1983-2009 is far smaller than what is detectable by our observational systems. Extended Data Fig. 4a and b shows ALL cloud trends without the subtraction of the 60° S-60° N average change. They exhibit patterns similar to those seen in Figs 1c and 2d.

LETTER RESEARCH
Another robust prediction by climate models is the rising height of high cloud tops at all latitudes 10,12,18,20 . Figure 3a displays ISCCP climatological zonal mean cloud amount within 7 cloud top pressure intervals. Only amounts of clouds with optical thickness greater than 3.6 are plotted to reduce uncertainty in cloud-top pressure retrievals. A local maximum in cloud amount occurs in the 180-310 hPa interval in the tropics, whereas clouds are typically no higher than 310 hPa in the mid-latitude storm tracks, following the latitudinal variation of tropopause height. Figure 3b, c shows that the ISCCP and PATMOS-x zonal mean cloud amounts increased in the 50-180 hPa interval and decreased in the 180-310 hPa interval during the period 1983-2009 in the tropics, consistent with a rise in the tops of the highest clouds. The increase in cloud amount in the 180-310 hPa interval at mid-latitudes similarly suggests a rise in the highest cloud tops. Figure 3d displays trends in zonal mean cloud amount during the period 1983-2009 from the ALL ensemble mean. The pattern of modelled cloud trends is highly correlated with the satellite record in the 50-180 hPa and 180-310 hPa intervals, suggesting that the observed rise in cloud top is at least partly due to external radiative forcing.
We expect less agreement below these levels because the ISCCP and PATMOS-x satellite retrievals cannot detect lower clouds underneath higher clouds, whereas the models report the exact cloud amount at each level. We note that the negative trends in cloud amount occurring in the 50-180 hPa interval at higher latitudes are relative to the 60° S-60° N average cloud change for that interval and do not correspond to an actual reduction in cloudiness. Extended Data Fig. 4c, for which the 60° S-60° N average change was not subtracted, shows that modelled cloud amount in the 50-180 hPa interval merely increases less at higher latitudes than at lower latitudes.
What specific factors are contributing to the observed cloud changes? We addressed this question by examining the additional CMIP5 simulations listed in Extended Data Table 1 with external radiative forcing only from greenhouse gases (GHG), only from anthropogenic aerosol (AA), only from ozone (OZ), and only from natural solar variations and volcanic aerosol (NAT). Extended Data Figs 5, 6 and 7 display the ensemble mean modelled cloud trends for GHG, AA, OZ and NAT simulations. The GHG and NAT simulations both produce modelled cloud trend patterns that are significantly correlated with the observed cloud trend pattern (Table 1 and Extended Data Fig. 3). This includes decreasing total cloud amount at mid-latitudes (GHG and NAT), increasing total cloud amount in the northern tropics (NAT), and increasing cloud amount in the 50-180 hPa interval in the tropics and in the 180-310 hPa interval at mid-latitudes (GHG and NAT). In contrast, the AA and OZ simulations do not produce cloud trends that globally resemble the observed cloud trends, as demonstrated by insufficiently positive correlation ( Table 1). The OZ simulations do exhibit reduced cloud amount at Southern Hemisphere mid-latitudes 21 .
Both the GHG and the NAT simulations experience increasing tropospheric temperature and decreasing stratospheric temperature from the 1980s and the 2000s. This is caused by increasing greenhouse gases in the former case and a recovery from the 1982 El Chichón and 1991 Pinatubo volcanic aerosol episodes in the latter case [22][23][24] . Tropospheric warming and stratospheric cooling promote an increase in the height of the highest cloud tops 25,26 , and together with global warming, promote an expansion of the tropical zone and a poleward shift of storm tracks 27,28 . Depleted stratospheric ozone is an additional factor promoting a poleward shift of the Southern Hemisphere storm track 21,29 .
The expansion of subtropical dry zones results in less reflection of solar radiation back to space. As cloud tops rise, their greenhouse effect becomes stronger. Both of these cloud changes have a warming effect on climate. Our results suggest that radiative forcing by a combination of anthropogenic greenhouse gases and volcanic aerosol has produced observed cloud changes during the past several decades that exert positive feedbacks on the climate system. We expect that increasing greenhouse gases will cause these cloud trends to continue in the future, unless offset by unpredictable large volcanic eruptions.
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

Satellite datasets.
ISCCP provides values of cloud amount in seven cloud top pressure intervals and six optical thickness intervals (that is, cloud amount for each of 42 'cloud types') from July 1983 to December 2009 (ref. 6). Total cloud amount is the sum over all intervals. Cloud-top pressure is most accurately identified when clouds are nearly opaque at thermal infrared wavelengths. This occurs when cloud optical thickness at visible wavelengths is greater than 3.6. Geostationary satellites are the primary contributors to the ISCCP cloud record. We downloaded ISCCP D1 data from the Atmospheric Science Data Center located at NASA Langley Research Center and applied a correction procedure to remove spurious variability associated with changes in satellite orbits, satellite calibration and ancillary data 5 . We note that the correction procedure removes any real global mean cloud variability, so all trends presented in this study are with respect to an unknown global mean trend, which could be zero. The present study uses only daytime observations (defined as solar zenith angle <78°) because visible radiances are required to retrieve cloud optical thickness. We found that trends in total cloud amount retrieved from day and night infrared radiances are very similar to trends in total cloud amount retrieved from daytime visual+infrared radiances.
PATMOS-x provides cloud amount in seven cloud-top pressure intervals and six optical thickness intervals starting in October 1981 (ref. 7). The present study uses data from January 1983 to December 2009 for consistency with ISCCP. The total cloud amount is the sum over all intervals. Polar-orbiting satellites are the only contributors to the PATMOS-x cloud record. We downloaded PATMOS-x Version 5 Level 3 'GEWEX' data from https://cimss.ssec.wisc.edu/patmosx/access. html and applied a correction procedure to remove spurious variability associated with changes in satellite orbits, satellite calibration, and ancillary data 5 . As for ISCCP, the correction procedure removes any real global mean cloud variability. For consistency over the entire PATMOS-x record, we use products retrieved only from the daytime pass of the 'afternoon' satellites.
The passive remote sensing techniques employed by ISCCP and PATMOS-x can have difficulty identifying the occurrence of optically thin cirrus overlying optically thick lower cloud and can underestimate the height of the cloud top when cloud particle density is sparse in the upper portion of an optically thick cloud 30,31 . ISCCP suffers more from remote sensing limitations than PATMOS-x since the latter dataset uses more wavelengths and thus has more information available to retrieve cloud properties 7 . The result is a downward bias in reported cloud-top height relative to that obtained from active remote-sensing techniques, for which only a short record is available. ISCCP correspondingly underestimates the cloud amount in the 50-180 hPa and 180-310 hPa pressure intervals compared to active remote sensing. Despite the bias, a real increase over time in the height of the highest cloud tops will nonetheless be reported by ISCCP as an increase in the amount of cloudiness in the higher elevation pressure interval and a decrease in the lower elevation pressure intervals. The bias may produce an underestimate in the magnitude of cloud trends since ISCCP climatological cloud amount is underestimated, but this does not undermine our analysis because we compare the relative spatial patterns of observed and modelled cloud change rather than the absolute magnitudes of observed and modelled cloud change.
Albedo is a useful parameter for our investigation because variability in cloud amount is by far the dominant contributor to variability in albedo outside icecovered regions. Variability in cloud optical thickness is a secondary contributor. ERBS albedo values are available for November 1984 through February 1990, but we use the January 1985-December 1989 climatology for simplicity 13 . We also use measurements from ERBS only because the other two satellites contributing to the Earth Radiation Budget Experiment (NOAA9 and NOAA10) were not available for the entire period. We note that ERBS was in a precessing orbit and sampled the entire diurnal cycle. CERES Energy Balanced And Filled (EBAF) Ed2.8 reflected solar radiation values are available from the morning satellite Terra starting in March 2000 and from the afternoon satellite Aqua starting in July 2002 (ref. 14). Since Terra and Aqua are in Sun-synchronous orbits CERES EBAF uses geostationary satellites to fill out the diurnal cycle of cloudiness not sampled by Terra and Aqua 15 . To ensure consistent sampling of the diurnal cycle and seasonal cycle, we constructed a climatology for July 2002-June 2014. We then divided reflected solar radiation by incoming solar radiation at the top of the atmosphere to calculate albedo. ERBS data were obtained from a CD-ROM provided by the Atmospheric Science Data Center located at NASA Langley Research Center, and CERES data were obtained from the NASA Langley Research Center CERES ordering tool at (http://ceres.larc.nasa.gov/). Although the datasets are individually well calibrated, there is no absolute calibration between ERBS and CERES. To bring them to a common reference point, we multiplied the ERBS albedo by a constant factor so that ERBS and CERES had the same climatological annual albedo averaged over 60° S-60° N. This means that CERES-ERBS differences are relative to an unknown global mean difference that could be zero. Version 4 of the MAC-LWP dataset for January 1988-December 2014 provides a useful complement to measurements of cloud amount and albedo 16 . The MAC-LWP dataset synthesizes passive microwave retrievals from 12 different sensors using the Remote Sensing Systems version 7 ocean algorithm 32 . Liquid water path is the spatially averaged vertically integrated amount of cloud liquid water within a satellite footprint. Cloud-free areas contribute a value of zero to the spatial average within the footprint. Liquid water path increases as clouds become more horizontally extensive (that is, as cloud amount increases). It also increases as clouds become vertically thicker and as cloud water concentration becomes larger. The dataset does not include contributions from cloud ice, and retrievals are available only over open ocean. To provide a similar basis for comparison to the ISCCP, PATMOS-x, and ERBS/CERES datasets, from which global mean variability was removed in the correction and adjustment process, we subtracted the 60° S-60° N average liquid water path from the value at each grid box for each month. This has little impact on the spatial distribution of trends. Total cloud amount is available from 107 realizations from 33 models, and cloud amount in each vertical layer is available from 76 realizations from 27 models (Extended Data Table 1). We calculated the ensemble mean as a simple average of all available realizations. Some models provided only one realization and other models provided up to ten realizations for the same external forcing. Natural internal variability across the simulations tends to cancel in the ensemble mean, leaving behind the radiatively forced component of cloud change. The ensemble mean has smaller trend amplitudes than any one realization or the observations due to this reduction of natural internal variability.
A smaller set of CMIP5 models provided additional simulations with external radiative forcing only from anthropogenic greenhouse gases (GHG), only from anthropogenic aerosol (AA), only from ozone (OZ), and only from natural solar variability and volcanic aerosol (NAT) (Extended Data Table 1). A few models included ozone variability in GHG simulations, but we excluded these from our analysis to avoid confusion about forcing factors. Total cloud amount is available from 44 realizations from 14 models for GHG, from 33 realizations from 11 models for AA, from 11 realizations from 3 models for OZ, and from 37 realizations from 10 models for NAT. Cloud amount in each vertical layer is available from 35 realizations from 12 models for GHG, 32 realizations from 11 models for AA, from 1 realization from 1 model for OZ, and 28 realizations from 8 models for NAT. We did not analyse cloud amount in vertical layers for OZ since only one realization was available.
For many of these models, GHG, AA, OZ, and NAT output was available only until 2005. To maximize the number of realizations, we chose the 1979-2005 interval to calculate trends for GHG and AA since this time period has the same length as the ISCCP and PATMOS-x records. The four-year shift for trend calculation should not matter for the GHG and AA simulations since greenhouse gas and aerosol emissions vary on multidecadal rather than interannual timescales. Timing matters more for OZ and NAT because the former includes stabilization of the ozone hole in the 2000s and the latter includes volcanic eruptions, so we chose only those models providing output for the full 1983-2009 time period. We note that the set of contributing models and numbers of realizations is not identical for the ALL, GHG, AA, OZ, and NAT simulations. We chose to use all available simulations from each forcing scenario because restricting our comparison to only those models and numbers of realizations in common would result in a much smaller sample size.
Most CMIP5 models provided multicentury simulations of preindustrial conditions with no anthropogenic or natural external radiative forcing as a control case (Extended Data Table 2). Cloud variability in these simulations results only from natural internal variability of the coupled ocean-atmosphere-land system, including El Niño/Southern Oscillation (ENSO). Total cloud amount is available for a total of 15,807 years from 27 models, and cloud amount in each vertical layer is available for a total of 7,311 years from 13 models.
Output from CMIP5 simulations was downloaded from the Earth System Grid Federation. To provide a similar basis for comparison to the ISCCP, PATMOS-x, and ERBS/CERES datasets, from which global mean variability was removed in the correction and adjustment process, we subtracted the 60° S-60° N ocean average cloud amount from the value at each ocean grid box for each month and subtracted the 60° S-60° N land average cloud amount from the value at each land grid box.
With the exception of cloud amount in the 50-180 hPa interval, this has little impact on the spatial distribution of trends. For zonally averaged cloud amount in each vertical layer, we subtracted the 60° S-60° N average.
Since the CMIP5 models do not routinely report amount of cloudiness in various optical thickness intervals, we could not limit our analysis of modelled clouds in Fig. 3 to only those with optical thickness greater than 3.6. Cloud amount from the satellite records also differs from standard cloud amount output from CMIP5 models in that the former do not detect clouds with optical thickness less than about 0.3, whereas the latter report the amount of all clouds, even very optically thin ones. Another difference is that CMIP5 models report actual cloud amount at each model layer, not cloud amount unobscured by higher clouds as do the ISCCP and PATMOS-x satellite datasets. For a closer comparison to observations, a few CMIP5 models incorporated the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) 33 . This software produced model cloud output according to how it would be detected through the limitations of satellite retrieval, most often ISCCP. Since it is computationally expensive, COSP cloud output is available only for 13 ALL realizations from seven models (Extended Data Table 1). We use standard model cloud output in order to have a larger sample size but obtain similar but noisier results if we use COSP cloud output from an ISCCP satellite simulator (Extended Data Fig. 8). Data analysis. ISCCP, PATMOS-x, and ERBS data on a 2.5° × 2.5° latitudelongitude grid and CERES and MAC-LWP data on a 1° × 1° grid were spatially averaged to a common 5° × 5° grid. Output from CMIP5 models on a variety of grid sizes were bilinearly interpolated to a 5° × 5° grid. We linearly interpolated cloud amount from model layers to a common vertical grid with regular 50 hPa spacing in pressure. To display modelled cloud trends in the figures, we linearly interpolated trends from the 50 hPa regular grid to the midpoints of the seven pressure intervals used by the satellite datasets. All calculations are performed on anomalies from the long-term mean. In all cases of spatial averaging and spatial correlation, grid box values were weighted by the cosine of the grid box centre latitude to account for the variation of grid box area with latitude. We restrict our analysis to latitudes equatorward of 60° because passive retrieval of cloud properties by satellite is difficult over bright and cold surfaces, and no visible retrievals can be made during polar night.
We use least-squares linear trends or the average difference between two time periods as convenient means of summarizing change over time. Two-sided P values for the trends are determined using a conventional Student's t-test with an effective sample size that takes temporal autocorrelation into account.
For simplicity of comparison with modelled total cloud amount trends in calculating correlations in Table 1, we averaged cloud/albedo changes from all cloud amount and albedo datasets together before comparing to the ensemble mean total cloud amount trends from the CMIP5 ALL, GHG, AA, OZ, NAT and preindustrial simulations. Since cloud amount and albedo have different physical units, we standardized the grid box changes before averaging. Specifically, we divided each ISCCP grid box cloud amount trend by the standard deviation of all ISCCP grid box cloud amount trends, each PATMOS-x grid box cloud trend by the standard deviation of all PATMOS-x grid box cloud amount trends, and each grid box albedo change by the standard deviation of all grid box albedo changes. For simplicity of comparison with modelled cloud amount trends in seven pressure intervals, we averaged cloud trends from ISCCP and PATMOS-x datasets together.
Our null hypothesis is that the observed cloud changes result purely from natural internal variability. If so, there should be no systematic relationship between the spatial pattern of cloud trends generated by natural internal variability and the spatial pattern of cloud trends generated by external radiative forcing. The former is represented by individual preindustrial simulations, each with different realizations of natural internal variability, and the latter is represented by the ensemble mean of simulations with external radiative forcing, in which natural variability has been largely averaged out. The suitability of this null hypothesis can be demonstrated by calculating the distribution of spatial correlation coefficients (Pearson's r) between the pattern of cloud trends from the ensemble mean of forced simulations (ALL, GHG, or NAT) and the pattern of cloud trends from time periods of similar length obtained from preindustrial control simulations. We build up each null distribution by calculating cloud trends and spatial correlation values during a rolling 27-year period (that is, years 1-27, years 2-28, years 3-29, and so on) through the preindustrial control simulation for each model. Extended Data Fig. 3 displays the frequency distributions of the calculated correlation values. There are 15,104 time periods for total cloud amount and 6,973 time periods for cloud amount in vertical layers. The mean and median correlation values of the null distributions are zero, as expected, and the total area under each frequency distribution is equal to unity.
Our alternative hypothesis is that external radiative forcing was a contributing factor in producing the observed cloud trends. If so, we expect a positive spatial correlation between the observed trend pattern and the trend pattern from the ensemble mean of simulations with external radiative forcing (values shown as vertical lines in Extended Data Fig. 3). The P value for a particular spatial correlation value r is simply the fraction of correlation values from the preindustrial control simulations with values more positive than r (that is, the fraction of area under the frequency distribution to the right of the vertical line). For simplicity, we calculate P values with respect to a null distribution built from spatial correlation values for cloud trends from single time periods. Another option is to build a null distribution from spatial correlation values for cloud trends from ensemble means of multiple, randomly-selected time periods. Our results are the same for either approach.
For corroboration of P values calculated with respect to cloud trend patterns from the preindustrial control simulations, we additionally computed one-sided P values using a conventional Student's t-test for statistical significance of Pearson's r. In this case, a critical parameter is the effective number of spatial degrees of freedom 34 . We determined that there are 51 spatial degrees of freedom between 60° S and 60° N for the observed total cloud amount in grid boxes. Simplistically considered, this corresponds to a set of boxes slightly larger than those outlined by the latitude-longitude grid lines in Fig. 1, if apportioned equally. However, remote teleconnections contribute to reduced spatial degrees of freedom in addition to local spatial coherence. Zonal means have substantially fewer degrees of freedom (8 for total cloud amount and 7 for cloud amount in the 50-180 hPa and 180-310 hPa pressure intervals). This corresponds to about 15° spacing in latitude. P values obtained from formal tests are in some cases substantially larger than those obtained from the preindustrial control simulations (Table 1 and Extended Data Fig. 3), suggesting that the actual number of effective spatial degrees of freedom may be larger than that indicated by the method we used 34 . ENSO-like variability. The dominant source of multiyear natural variability in the climate system is the El Niño/Southern Oscillation (ENSO) phenomenon. Variability occurring at interdecadal time scales, especially over the Pacific basin, exhibits a pattern similar to that of ENSO 35 . We investigated whether the observed cloud changes are a manifestation of ENSO-like variability by calculating the correlation of the spatial pattern of cloud trends with the spatial pattern of the difference between observed La Niña composite cloud anomalies and El Niño composite cloud anomalies (figure not shown due to space limitations). The correlation between the observed La Niña-El Niño pattern and the observed trend pattern is only 0.13, and the spatial correlation between the observed La Niña-El Niño pattern and the ensemble mean ALL trend pattern is only 0.14. Considering that the spatial correlation between the observed trend pattern and the ensemble mean ALL trend pattern is 0.39 (Table 1), we think ENSO-like variability cannot be a major contributor to the global pattern of cloud change between the 1980s and the 2000s. Code availability. All plots and calculations were produced with custom code using the NCAR Command Language (NCL) Version 6.3.0 (http://dx.doi. org/10.5065/D6WD3XH5). This code can be obtained by contacting the corresponding author.