Cloud responses to climate variability over the extratropical oceans as observed by MISR and MODIS

Linear temporal trends in cloud fraction over the extratropical oceans, observed by NASA's Multi-angle Imaging SpectroRadiometer (MISR) during the period from 2000 to 2013, are examined in the context of coincident European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data using a maximum covariance analysis. Changes in specific cloud types defined with respect to cloud-top height and cloud optical depth are related to trends in reanalysis variables. A pattern of reduced high-altitude optically thick cloud and increased low-altitude cloud of moderate optical depth is found to be associated with increased temperatures, geopotential heights, and anti-cyclonic flow over the extratropical oceans. These and other trends in cloud occurrence are shown to be correlated with changes in the El Niño–Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), the North Pacific index (NPI), and the Southern Annular Mode (SAM).

1 Introduction

Clouds play a fundamental role in Earth's climate due to their effect on the planet's radiative budget. However, cloud responses to climate change are poorly understood, and cloud–climate interaction is presently one of the largest sources of uncertainty in climate models (Caldwell et al., 2016; Boucher et al., 2013; Bony et al., 2006). Several changes in midlatitude cloud are expected under global warming with varying degrees of certainty: including poleward shifts in the storm tracks, rising melting level, rising high cloud tops, and reduced low cloud (Boucher et al., 2013). Understanding changes in midlatitude and Southern Ocean cloud in particular is important, because these clouds have a large radiative impact, influence atmospheric dynamics (Kay et al., 2016; Hwang and Frierson, 2013), and are not adequately captured by climate models (Trenberth and Fasullo, 2010; Bodas-Salcedo et al., 2014; Kay et al., 2012). Several studies have observed midlatitude cloud responses to extratropical synoptic variability, for instance, changes in cloud cover associated with the North Atlantic Oscillation (NAO) or Southern Annular Mode (SAM) (e.g. Li et al., 2014a, b; Li and Thompson, 2016; Ceppi and Hartmann, 2015; Gordon and Norris, 2010; Gordon et al., 2005; Tselioudis et al., 2000). Many studies of cloud variability (including several of those cited above) are based on analysis of International Satellite Cloud Climatology Project (ISCCP) data sets. ISCCP is a multi-instrument and multi-satellite product that combines observations from polar orbiting and geostationary weather satellites to determine cloud amount and categorizes clouds by their cloud-top pressure and optical depth (Rossow and Schiffer, 1999). For example, Bender et al. (2011) use meridional maxima in the ISCCP total cloud fraction as a proxy for the midlatitude storm track latitude, and identify a 25-year poleward trend in the storm track location. Building on early studies, Norris et al. (2016) examine trends in the ISCCP and Pathfinder Atmospheres–Extended data sets after applying several empirical corrections (Norris and Evan, 2015). They also show a poleward trend in the storm track location between the 1980s and 2000s and identify an increase in global cloud-top height. These findings are consistent with Hadley cell expansion and increased tropopause height predicted by climate models under increased CO₂ (Boucher et al., 2013).

While these studies are insightful, the ISCCP data suffer from difficulties related to changes in instrumentation, orbital configurations, viewing angle, and individual instrument calibration drifts that make it difficult to use these data to identify subtle changes in cloud fraction, especially for particular cloud types. Authors of these previous studies recognized and used various approaches to account for limitations of the data (Evan et al., 2007; Bender et al., 2011; Norris and Evan, 2015). Only recently have higher-quality cloud data sets from NASA's Earth Observing System (EOS) reached sufficient duration to begin to observe trends on intra-decadal timescales. Marchand (2013) found that both the Multi-angle Imaging SpectroRadiometer (MISR) and the MODerate resolution Imaging Spectroradiometer (MODIS) data sets contain linear temporal trends in cloud fraction (for specific cloud types) in various regions between 2000 and 2010. These trends include (a) a decrease in mid-level clouds over the Southern Ocean (30–50^∘ S and strongest over the South Pacific), (b) a decrease in high-level clouds and increase in low-level clouds over the North Pacific, and (c) an increase in cloud of moderate optical depth (1.3 < OD < 23) and an approximately compensatory reduction in cloud of high optical depth (OD > 23) over the extratropical oceans of both hemispheres.

To better understand the MISR and MODIS extratropical cloud trends, we explore linkages between the MISR extratropical cloud fraction trends and trends in ERA-Interim reanalysis variables. We do this using a maximum covariance analysis (MCA) to identify spatially coincident trends in cloud occurrence and several reanalysis variables. Furthermore, we examine correlations with known modes of climatic variability including the El Niño–Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), and the SAM.

2 Data

2.1 Cloud observations

We use the Multi-angle Imaging SpectroRadiometer (MISR) L3 1^∘ gridded CTH-OD Version 6 data set, from 2000 to 2013 (Marchand et al., 2010). MISR is a polar orbiting instrument aboard EOS Terra. MISR images the Earth using an array of nine push-broom style cameras with four spectral bands (red, green, blue, and near IR) oriented at different viewing angles along its orbit, and uses multi-view geometry to estimate cloud-top heights (and several other products) (Diner et al., 1998). In this data set, total cloud cover (or total cloud fraction) in each spatial grid cell is divided into a cloud-top height vs. cloud optical depth (CTH-OD) joint histogram. Cloud longwave radiative effects are directly related to cloud-top temperature which is related to cloud-top height, whereas cloud shortwave radiative effects are directly related to cloud optical depth. As such, the CTH-OD joint histograms are convenient for studying cloud–climate interaction because they categorize clouds based on physical properties which relate to their expected radiative effects. Here, “cloud fraction” is defined as the fractional occurrence of cloud (for any histogram bin or combination of bins) at a given location in a month relative to the total number of observations (pixels) observed by MISR. MISR observes most of the Earth every 8 days, and has a resolution (at nadir) of about 275 m from which cloud occurrence (detection), cloud-top height and cloud optical depth are determined over ice-free ocean on a 1.1 km grid. A detailed description of the MISR CTH and OD retrievals can be found in Marchand et al. (2010). Each L3 cell in the monthly 1^∘ grid typically has several thousand observations associated with it, although many of these will have been taken concurrently (with a high degree of correlation). At 1^∘ resolution, in the extratropics, the MISR cloud occurrence data have a zonal autocorrelation length of less than 5^∘, and we perform our analysis using 5^∘ spatial averages of the MISR data (such that each of our grid cells is effectively an independent sample).

Trends computed from MODerate resolution Imaging Spectroradiometer (MODIS) Collection 6, 1^∘ gridded cloud-top pressure and optical-depth histograms (Hubanks et al., 2015) are briefly compared with MISR trends. MODIS is a 36-band polar orbiting radiometer aboard both EOS Aqua and Terra, with equatorial crossing times of about 10:30 and 13:30 LT, respectively. MODIS provides cloud fraction joint histograms that are similar in structure to the MISR joint histograms, but are retrieved using different algorithms, which are described by Platnick et al. (2015) and compared with MISR and ISCCP retrievals in Marchand et al. (2010). In MODIS Collection 6 processing (Platnick et al., 2017), pixels that are determined to be only partly cloudy or on the edge of a cloud (meaning cloudy pixels that border a clear pixel) are stored in a separate histogram from other cloudy pixels, which are nominally fully cloudy. Inclusion of the partly cloudy and edge pixels does not have a substantial effect on the MODIS cloud fraction trends. Because the quality of partly cloudy and edge pixel retrievals is suspect (and for consistency with Marchand (2013) who used the MODIS Collection 5 cloud product that does not include partly cloudy or edge pixels), we show results without the partly cloudy or edge pixels in later figures.

2.2 Reanalysis

The ERA-Interim reanalysis (Dee et al., 2011) is maintained by the European Centre for Medium-Range Weather Forecasts (ECMWF). The data used here span the same period as the MISR data (March 2000–March 2013) and are the ERA-Interim “monthly means of daily means.” The MISR cloud fraction data are compared to several reanalysis variables: temperature (T), geopotential (Z), specific humidity (q), vertical velocity (ω), divergence (∇), absolute vorticity (ζ), sea surface temperature (SST), and sea level pressure (SLP). The data used are defined at 20 pressure levels: 1000–50 hPa by increments of 50. Because MISR retrieves physical cloud-top height, as opposed to cloud-top pressure, and because the ECMWF data are defined at pressure levels, the hypsometric equation was used to derive altitudes (z) for the monthly ECMWF data, and the reanalysis data were then linearly interpolated onto MISR's cloud-top height grid.

It should be noted here that potentially spurious features exist in the reanalysis specific humidity trends. The ERA-Interim data indicate a ubiquitous increase in low-level tropospheric specific humidity nearly everywhere, but particularly in the tropics, which is in the order of 0.1–0.3 g kg⁻¹ decade⁻¹, and occurs primarily below 750 hPa. We performed a brief comparison to specific humidity data from the Modern-Era Retrospective analysis for Research and Applications (MERRA) and found that while it did indicate an increase in low-level specific humidity of similar magnitude in some of the Northern Hemisphere, it did not corroborate the pervasive trend in the ECMWF data. Dessler and Davis (2010) provide a more comprehensive intercomparison of specific humidity trends in different reanalysis data sets, albeit not for the time period analyzed here, and conclude that while most show recent increases in specific humidity (which is expected under global warming), there is large disagreement over the magnitude and spatial distribution. In any case, we show the ECMWF specific humidity data, but caution the reader that these data may not be as robust or reliable as other fields.

2.3 Climate indices

Several of the patterns found via the MCA (Sects. 4–6) resemble well-documented modes of climate variability. In Sect. 4, we compare the monthly time series of the MCA modes to monthly indices for various modes of climate variability and Northern Hemisphere teleconnection patterns maintained by the NOAA Climate Prediction Center (CPC). The indices used are as follows:

• The Niño region 1 + 2, 3, 3.4, and 4 indices (hereafter referred to as “Niño 3.4,” for instance). These indices are based on spatial averages of SST anomalies in various regions in the tropical Pacific Ocean: 1 + 2 refers to a latitude of 0–10^∘ S and longitude of 90–80^∘ W; 3 refers to a latitude of 5^∘ N–5^∘ S and longitude of 150–90^∘ W; 3.4 refers to a latitude of 5^∘ N–5^∘ S and longitude of 170–120^∘ W; and 4 refers to a latitude of 5^∘ N–5^∘ S and longitude of 160^∘ E–150^∘ W.

• The PDO index, which is the first mode of an empirical orthogonal function (EOF) decomposition of SST north of 20^∘ N in the Pacific Ocean (Mantua et al., 1997).

• The Antarctic Oscillation index or “Southern Annular Mode” (SAM), which is defined as the first mode of an EOF analysis of 700 hPa geopotential height south of 20^∘ S 1979–2000 (Thompson and Wallace, 2000).

• The Pacific–North American (PNA) mode index. The index is defined by the projection of the PNA loading pattern onto the daily 500 hPa height anomalies over the entire Northern Hemisphere. The PNA loading pattern is derived by a rotated principal component analysis of 500 hPa heights north of 0^∘ between 1950 and 2000 as described in Barnston and Livezey (1987).

• The North Pacific index (NPI). The NPI is a standardized mean of sea level pressure between 30–65^∘ N and 160–220^∘ E (Trenberth and Hurrell, 1994), and was obtained from the University Corporation for Atmospheric Research (UCAR) website (link provided in the data availability section).

Several other similarly defined teleconnection patterns were analyzed, but are less relevant to our results: the East Pacific/North Pacific, Scandinavian, Tropical/Northern Hemisphere, East Atlantic, Pacific Transition, Polar/Eurasia, and West Pacific indices. All indices but the NPI were obtained from the NOAA CPC website (see data availability section). These indices are dimensionless except for the “Niño” indices, which represent temperature anomalies, although we standardized each index prior to computing any statistics by subtracting the mean and dividing by the standard deviation. Correlation coefficients computed using these indices (shown in Fig. 4) are computed after first detrending the indices, the removed trends are shown in Fig. 4b.

Figure 1 Cloud fraction trends (2000–2013), in percentage per decade, computed for each bin of MISR cloud-top height vs. optical-depth joint histograms (a)–(d). Trends are spatially averaged over the four extratropical ocean basins listed above each panel (and these regions are depicted in Figs. 2 and 3). Bold bordered bins denote trends that are significant at the 95 % level. The middle set of panels (e)–(h) show spatially averaged trends in cloud fraction with respect to the optical-depth bins only (summing over all cloud-top height bins), with whiskers denoting 95 % confidence intervals. The bottom panels (i)–(l) show similar cloud optical-depth trends for MODIS Aqua (blue) and Terra (red). The MODIS CTH-OD product uses slightly different optical-depth bins and two of the MODIS high optical-depth bins have been summed over to create one OD > 60 bin for easier comparison with the MISR trends. MISR “no retrieval” bins, denoted “NR,” include cases where either the cloud-top height or optical-depth retrieval failed. MISR observed a significant reduction in cloud fraction in one or both of the τ>23 bins in each ocean basin, and three (all but the South Atlantic) show a significant increase in low-level cloud of moderate optical thickness. MODIS Aqua is largely consistent with MISR (it shows increases in clouds of moderate optical thickness in the same three basins, and differs mostly in that it shows no statistically significant decrease in clouds with an optical thickness greater than 23 in the South Pacific). MODIS Terra is less consistent with MISR (and MODIS Aqua) and shows smaller decreases or slight increases in clouds with an optical thickness greater than 23, likely due to calibration issues (see Sect. 2 of text for discussion).

Figure 2 Trends between 2000 and 2013 in ERA-Interim 500 hPa geopotential height, temperature, and absolute vorticity in each of the extratropical study regions: North Atlantic (25–65^∘ N, 280–360^∘ E), North Pacific (25–65^∘ N, 120–240^∘ E), South Atlantic (25–65^∘ S, 280–360^∘ E), and South Pacific (25–65^∘ S, 120–240^∘ E). Hatching denotes trends that pass a 95 % confidence test. The North Pacific region shows significantly increased heights (d), temperature (e), and anti-cyclonic flow (f) (negative trends in absolute vorticity in the Northern Hemisphere) primarily in the storm track regions and extending into the subtropics. All of these features are commonly associated with extratropical high-pressure systems. Both of the Southern Hemisphere study regions show similar increases in heights (g, j), temperature (h, k), and anti-cyclonic flow (i, j) in the storm track regions as well; however, only small areas in the center of the locations with positive trends pass the significance test. The North Atlantic study region (a–c) shows no trends that are significant at the 95 % level.

3 Trends in cloud fraction

Linear temporal trends were computed on the MISR cloud fraction data for the 2000–2013 period. Figure 1a–d show the trends computed in each of four extratropical ocean basins for each bin in MISR's CTH-OD joint histograms. The four regions studied are the North Atlantic (25–65^∘ N, 280–360^∘ E), the North Pacific (25–65^∘ N, 120–240^∘ E), the South Atlantic (25–65^∘ S, 280–360^∘ E), and the South Pacific (25–65^∘ S, 120–240^∘ E), (these geographic regions are shown in Figs. 2 and 3, which display trends in ERA-Interim data and results of the MCA that is discussed in Sect. 4). The cloud fraction data were spatially averaged over each ocean basin prior to computing trends, and the composited seasonal cycle was removed. Figure 1e–h show the MISR cloud fraction trends associated with each optical-depth category, where cloud fraction is summed with respect to cloud-top height prior to computing trends. Finally, Fig. 1i–l show cloud fraction trends for each MODIS optical-depth category for both MODIS Aqua and Terra. Note that the MODIS cloud occurrence histograms do not include a “no retrieval” (NR) category. The MODIS data also include separate bins for optical depths between 60–100 and 100–150 (the high optical-depth bin is new for MODIS Collection 6); however, in Fig. 1, these two optical-depth bins have been summed to create a single bin representing all optical depths greater than 60. This step is taken to make the a comparison of the two data sets easier. Bold bordered bins in the joint histograms in Fig. 1a–d indicate that the cloud fraction trend in that bin exceeds a 95 % confidence test, whereas the bars in Fig. 1e–l indicate the 95 % confidence interval. The confidence intervals were computed using a windowed bootstrapping technique described in Wilks (2006), which was also used to assign confidence to the cloud fraction trends computed in Marchand (2013). This technique involves randomly resampling (with replacement) each bin's cloud-fraction time series in 12-month chunks 1000 times and computing trends for each of the resampled time series. The trend associated with the original time series is said to be significant at the 95 % level if it exceeds the 25th most positive or 25th most negative resampled trend. Bins that account for less than 0.1 % of the total cloud fraction are not considered. Figure 1 shows the same dominant pattern of changing extratropical cloud fraction identified in Marchand (2013), but here it is based on 3 additional years of MISR data. The trends in Fig. 1 do not have a strong seasonal dependence (not shown). The results in three of the basins (the North Pacific, South Atlantic, and South Pacific) are characterized by increasing cloud of moderate optical depth, particularly at low altitudes, and decreasing cloud of high optical depth at most levels (with the largest decreasing trends at high levels). The South Atlantic differs in that it does not have a clear increase in low-level clouds, only a reduction in optically thick clouds. The South Pacific also differs from the other basins, featuring a stronger reduction in mid-level cloud (between 2 and 5 km).

There are no significant changes in the occurrence of MISR failed optical-depth retrievals (the NR column), but there is a slight reduction in the number of failed cloud-top height retrievals (the NR row) in the North Pacific and North Atlantic. Failed cloud-top height retrievals most often occur in multilayer cloud conditions (where a low cloud that is visible in the MISR nadir view cannot be located in an off-nadir view due to a visibly opaque or semitransparent higher altitude cloud), which suggest that a small portion of the observed increase in low-level clouds is due to a reduction in higher (semitransparent) cloud.

We note here that an analysis of the MISR calibration suggests that MISR near-IR radiances (used to obtain MISR optical depths) likely have a small downward drift of about 0.9 % to 1.5 % per 10 years (Bruegge et al., 2014; Corbett and Loeb, 2015; Limbacher and Kahn, 2017). This calibration drift can be expected to reduce the retrieved optical depth, decreasing the occurrence of clouds with large optical depths in the CTH-OD product, and increasing the occurrences of clouds with moderate optical depths. Such a calibration drift will not change the cloud-top altitude but will cause clouds at a given altitude to shift toward lower optical depths at that same altitude. Evidence for such a calibration-driven change can be seen to some degree in Fig. 1, where in the South Pacific (Fig. 1d) and South Atlantic (Fig. 1c) between 5 and 7 km there is a strong increase of cloud with optical depths between 9.4 and 23, and strong decrease at this same altitude for optical depths greater than 60. However, Fig. 1 suggests that much of the reduction in optically thick cloud occurs at high altitudes, whereas the increase in clouds with moderate optical depths occurs for low-level clouds. Limited testing, where the CTH-OD data set was reprocessed for 1-month with the observed radiances reduced artificially by 2 %, suggests a reduction in the occurrence of cloud with large optical depths may well explain 50 % to 75 % of the MISR trend depending on the region. Plans are underway to reprocess the entire MISR mission, starting from the (level 1) calibrated imagery and eventually including all higher-level data sets (level 2 swath and level 3 global data sets), which also incorporates the CTH-OD product. This reprocessing will include a correction for this calibration drift among other issues. For the present, an important caveat is that the strength of the optical-depth trend in the MISR data set (and associated statistical confidence) is likely being overestimated, and this may explain why the trends in the MISR data set are larger than MODIS. However, the significant changes in cloud-top height observed in three of the study regions and the reduced mid-level cloud (between about 2 or 2.5 and 4 km) observed in the South Pacific could not have been caused by this calibration drift.

In Fig. 1i–l, trends in MODIS cloud fraction bins only partially corroborate those in the MISR data set. MODIS Aqua identifies a reduction in optically thick cloud in all of the regions studied except the South Pacific, whereas MODIS Terra shows a reduction in only the North and South Atlantic. MODIS Aqua shows an increase in cloud of moderate optical depth in all regions except the South Atlantic (in agreement with MISR), whereas MODIS Terra shows little or no change in these bins. As with MISR, there is evidence for drifts in the MODIS calibration (Lyapustin et al., 2014). Corbett et al. (2015) compares both MISR and MODIS Terra level 1 radiances to collocated CERES outgoing shortwave radiation observations and finds that while MISR red, green, and near IR bands have darkened relative to CERES, MODIS Terra red and near IR bands have brightened, which likely explains much of the discrepancy between MODIS Terra and MISR cloud optical thickness trends. Taken in combination, the MISR and MODIS data suggest that there has been a reduction in cloud optical thickness during the period studied, at in least in the North and South Atlantic Ocean and likely the North Pacific as well.

As another note, we add that there is a known issue with the MODIS Terra cloud mask over ocean. MODIS Terra's 8.6 µm channel has undergone warming since around 2010 that has not been corrected via onboard calibration, and this warming has caused a number of clear pixels to be flagged as cloudy in the MODIS cloud product. This problem primarily affects the low cloud retrieval fraction in tropical and subtropical regions with low average total cloud fraction and does not appear to have a substantial impact on the MODIS Terra extratropical cloud trends shown in Fig. 1. This trend is being corrected in MODIS Collection 6.1.

On monthly timescales cloud occurrence is heavily influenced by synoptic conditions, and it seems likely that a large portion of the observed cloud fraction trends are related to trends in synoptic variables. In Fig. 2, we show the spatial distribution of trends in several of the reanalysis variables discussed in Sect. 2 (500 hPa geopotential height, temperature, and absolute vorticity), hatching denotes trends that are significant with 95 % confidence (using 5×5^∘ bins). These trends have been computed after first deseasonalizing the data using compositing, and confidence intervals are determined using the windowed bootstrapping technique discussed above. There is a notable increase in both the 500 hPa temperature and geopotential height in the center of the North Pacific region (Fig. 2d–e). This is accompanied by increased anti-cyclonic flow (Fig. 2f) (a negative trend in absolute vorticity in the Northern Hemisphere), which is expected given the temperature and pressure changes, although the trends in cyclonicity often do not pass the significance test, which is perhaps due to the large variability in this field. Both the 500 hPa temperature and absolute vorticity fields show increasing trends in both Southern Ocean storm track regions (around 40–50^∘ S), with only small areas that pass the confidence test which coincide with the largest positive trends (Fig. 2h, i, k, l). These are accompanied by an increase in geopotential heights that does not pass the 95 % significance test, but is meteorologically consistent with the corresponding changes in temperature and vorticity (Fig. 2g, j). In the South Pacific this primarily occurs to the south and the east of New Zealand. In the South Atlantic there is a meridional dipole with the strongest positive changes in these fields in the poleward part of the region. Notably, the regions in the Southern Ocean that show increased geopotential heights, temperature, and anti-cyclonic flow correspond well with the loading pattern of the Southern Annular Mode which has undergone a positive trend during the study period and will be discussed in more detail later (Sect. 5.3 and 5.4). In the North Atlantic, while there is a weak increase in anti-cyclonic motion (c), as well as 500 hPa geopotential height (a) and temperature (b), along the storm track (which is similar in sign to the North Pacific) the trends are generally smaller than in the other three basins, with essentially no significant trends.

Figure 3 Results of the maximum covariance analysis between the MISR cloud fraction and the ECMWF reanalysis variable trends (2000–2013). The MCA modes are computed by applying a singular value decomposition to the spatial covariance matrix between MISR cloud fraction joint histogram trends and trends in vertical profiles of six ERA-Interim variables and sea surface temperature, as described in Appendix A. The middle column shows the pattern of changing cloud fraction in CTH-OD space (${\mathbf{U}}_{\mathbf{M}}^{*}$ in Appendix A) and the right column shows the associated trends in the vertical profiles of the reanalysis variables (${\mathbf{U}}_{\mathbf{E}}^{*}$). The left column shows the spatial loading patterns obtained by projection of the ECMWF modes on the far right onto the original trends. The spatial patterns are dimensionless, and the scales associated with the CTH-OD trend patterns are shown in the bottom right color bar, whereas the scales for the ERA-Interim profile trend patterns are shown in the legend near the top. The trend (due to each mode) at each latitude and longitude grid point is obtained by multiplication of the joint histogram or vertical profiles with the dimensionless spatial pattern on the left. Listed with each mode (in parentheses under the plot of the mode's spatial loading pattern) is the fraction of the covariance explained by that mode, followed by the fraction of variance in the ECMWF and MISR data sets, respectively. A 5^∘ latitude/longitude grid was used.

In summary, in the ocean basins studied, ERA-Interim reanalysis shows increased temperature, geopotential heights, and anti-cyclonic flow in the storm track latitudes, although more poleward in the South Atlantic. The largest and most robust changes are in the central North Pacific Ocean and weakest in the North Atlantic. Except perhaps in the North Pacific, most of the trend (linear fit to the change) in the reanalysis variables is not statistically significant; however, we stress that this does not mean that there has been no change in meteorology or that these changes have no impact on the clouds. To the degree the reanalysis data is accurate there has been an increase in temperature, geopotential heights, and anti-cyclonic flow in some portion of each of these basins over the period examined. Rather, the lack of statistical significance means that relative to the annual variability, the change is small and could be a result of annual variability rather than a true trend in the mean temperature, geopotential heights, and cyclonicity with time. In the following section we use a maximum covariance analysis to identify linkages between trends in the reanalysis variables and the MISR cloud fraction data. Finally, we show that the changes discussed above are consistent with recent trends in the NPI and the SAM and explore correlations with other climate indices.

4 Maximum covariance analysis

We wish to identify relationships between two high-dimensional data sets. The MISR data set is a function of time, latitude, longitude, cloud-top height, and cloud optical depth, whereas the ERA-Interim variables described in Sect. 2 vary in time, height, latitude, and longitude. In particular, we wish to identify which trends in cloud occurrence observed by MISR are related to trends in meteorology. To do this, we identify patterns of trending cloud occurrence, in the context of MISR cloud occurrence histograms, which tend to be co-located with patterns of trending meteorological variables using maximum covariance analysis (MCA), the results of which are shown in Fig. 3. The primary advantage of the MCA is that it allows for the simultaneous analysis of many different cloud and meteorological fields, and is thus well suited for an exploration of relationships between changes in clouds and changes in meteorology. It involves first computing a covariance matrix between the cloud and meteorological data, which is calculation of the covariance between each cloud-type and meteorological-variable pair, and then applying a singular value decomposition to the covariance matrix. The SVD identifies patterns (or modes) that describe the greatest covariance between the two data sets. This process avoids making any a priori assumptions about the most relevant meteorological fields or cloud types to analyze and can identify regionally specific interactions between clouds and meteorology. MCA falls into the family of other linear decomposition techniques such as SVD, principal component analysis, empirical orthogonal function analysis, and canonical correlation analysis (Bretherton et al., 1992; Hannachi et al., 2007; Von Storch and Zwiers, 2001). Typically, in climate science, SVD (often interchangeably referred to as empirical orthogonal function analysis and principal component analysis) is applied to a single time-varying field to identify important spatial patterns that covary in time. This involves decomposition of the field's time-covariance matrix. SVD was used in the atmospheric sciences as early as Lorenz (1956) and Kutzbach (1967), although it became significantly more popular in the 1980–1990s when it was used to identify a number of well-known modes of climate variability, such as the North Atlantic Oscillation (Barnston and Livezey, 1987; Wallace and Gutzler, 1981), the PDO (Mantua et al., 1997), and the Antarctic Oscillation (Thompson and Wallace, 2000). MCA involves applying SVD to a cross-covariance matrix computed between two data sets. It is similar to canonical correlation analysis which involves decomposition of a cross-correlation matrix. For instance, Prohaska (1976) applies SVD to the temporal correlation matrix computed between monthly mean sea level pressure and temperature fields. MCA has also been historically used to understand climate data, in particular, since about 2000, it has been heavily utilized to identify interactions between sea surface temperature and atmospheric variables (Czaja and Frankignoul, 1999; Liu et al., 2006; Frankignoul et al., 2011, and references therein). It has also been used in other areas of the earth sciences, for instance in sea-ice modeling (Dirkson et al., 2015) and the study of the structure of the Madden–Julian Oscillation (Adames and Wallace, 2014).

A complete mathematical description of the formulation used here, as well as a discussion of the significance of MCA modes, is given in Appendix A, and results are shown in Fig. 3. In short, a covariance matrix is computed that represents spatial covariance between trends in the MISR and trends in the ERA-Interim variables. A singular value decomposition (SVD) is applied to the covariance matrix. One set of singular vectors identified by the SVD then represents patterns (or modes) of cloud trends in CTH-OD space (shown in the middle column of Fig. 3), and the other set represents corresponding patterns in trends of the vertical profiles of the reanalysis variables (shown in the right column of Fig. 3). These patterns (or modes as we will call them in later sections) can be projected onto the original trend data to see corresponding spatial patterns (the left column of Fig. 3). In the left column of Fig. 3, bright red colors represent regions where the cloud occurrence and reanalysis trends shown in the panels to the right have occurred. The same is true for the blue colors except that in these regions the trends are of the opposite sign. For example, in North Atlantic mode 1 (the first row in Fig. 3) the red region around 45^∘ N in the leftmost panel indicates that there has been an increase in low cloud of moderate optical depth and a reduction in high thick cloud in the North Atlantic, and at the same location there has been increased geopotential heights, temperature, and anti-cyclonic motion at most levels. A limitation of the MCA is that a single MCA mode can only represent a single vertical structure (Fig. 3). Only when viewed holistically, combining multiple MCA modes, can the MCA account for variations in the vertical structure with latitude and longitude. In the discussion below, we chose to retain the first three modes in the case of the North Atlantic and Pacific Ocean basins and only two modes in the case of the Southern Hemisphere ocean basins. This choice is partly qualitative. As discussed below, the third MCA modes in the Northern Hemisphere each show very organized spatial structure and represent similar patterns of changing low cloud and low-level humidity. This was not the case in the Southern Hemisphere. In particular, the third MCA mode in the South Atlantic explains only 1 % of the spatial variance in the cloud fraction trends for that region. In Sect. 5, we discuss the MCA modes and their geographic patterns in detail.

Table 1 The fraction of the basin-averaged cloud occurrence trends in Fig. 1 that can be explained by each MCA mode (Fig. 3). Calculation is explained in Appendix A. Because most of the observed trends (Fig. 1) are changes in low cloud (CTH < 2.5 km) and in thick cloud (OD > 23), separate estimates are given for these categories. Notably, some MCA modes account for very little of the basin-averaged trend. This is due to spatial cancellation in the basin average. In NP1 for instance there are increases in low cloud in the eastern side of the basin that are largely balanced by decreases in low cloud on the western side (see dipole loading pattern Fig. 3 NP1).

Because the spatial loading patterns shown in Fig. 3 are derived by projection of the MCA modes onto the cloud and reanalysis trends, there is no requirement that summing across multiple modes will reproduce the entirety of the observed average trend in the basin (Fig. 1). However, we can calculate the amount of the basin-averaged cloud trend that can be explained by the cloud–meteorology relationship represented by each MCA mode. This is done by projecting the cloud trend patterns (middle column of Fig. 3) onto the spatial loading patterns derived from the reanalysis trends (left column in Fig. 3), taking the spatial mean, and subtracting the result from the observed spatially averaged cloud trends (see Appendix A). The results of this process are given in Table 1, along with the amount of the low cloud (cloud-top height < 2.5 km) and thick cloud (cloud optical thickness > 23) trend explained, and are discussed more in Sect. 5. The percentages in Table 1 are not necessarily cumulative across multiple modes, but do provide some idea of the ability of individual modes to explain the observed cloud occurrence trends in the basin mean. A useful feature of the MCA is that some modes highlight trends that mostly cancel in the spatial mean but represent important cloud–climate interaction. For instance, we will show later that the first mode in the North Pacific captures the influence of the PDO on North Pacific cloud, but this mode accounts for almost none of the spatially averaged (basin mean) trend because changes in the western part of the basin approximately cancel out changes in the east.

Unlike an EOF analysis, the original monthly data sets cannot be fully recovered from the MCA modes alone (information is lost when trends are computed and because the SVD is performed on the spatial covariance matrix between the two data sets). However, monthly time series associated with each mode can be derived by linearly projecting the MCA mode patterns shown in Fig. 3 onto the monthly cloud fraction anomaly data. As noted earlier, and discussed further in the next section, many, but not all, of the MCA modes resemble well-documented modes of climate variability. In many respects this is not surprising as midlatitude cloud is intimately linked to the synoptic meteorology; therefore, one might expect that changes in synoptic conditions captured by the CPC indices will be correlated with changes in cloud cover. We investigate the relationships between the MCA modes and various climate indices (introduced in Sect. 2) using a time-correlation analysis. In this analysis, both MCA-derived time series and the climate indices are first detrended before Pearson correlation coefficients are computed. The correlation values (which will be discussed further in the next section) are tabulated in Fig. 4, along with information about any linear temporal trends in the associated CPC index time series, which are shown in Fig. 4b. Many of the signals contain some trend, and the correlation coefficients are computed on detrended data to avoid spurious correlations which might occur just by virtue of both signals containing a trend. For instance, the PDO might show correlation with the South Atlantic modes just because both signals happened to either increase or decrease during the period studied, even if there is no physical connection. Instead we compute these correlations after first detrending, and show the strength of the trend that was removed in Fig. 4b. In effect, our correlation coefficients show the strength of the relationship based on covariability at monthly timescales. In Fig. 5 we show some specific examples of the detrended MCA-derived and CPC index time series for the North and South Pacific regions, to provide better intuition about the correlation coefficients in Fig. 4. These time series and correlation coefficients are also discussed in more detail in the next section in the context of individual MCA modes.

Figure 4 Correlation matrix between time series computed for each MCA mode shown in Fig. 3 and each of the CPC climate indices (a). The numerical values in (a) are correlation coefficients with the decimal point omitted, bold values indicate significance (p<0.01), coloring also denotes the sign and magnitude of the correlation. Only CPC indices with significant correlations with at least one MCA mode are included. The climate indices are sorted such that those that describe the most variance averaged across all MCA modes appear near the top. Panel (b) shows normalized trends in the CPC indices between 2000 and 2013, along with their magnitude expressed in standard deviations, showing for example, that the Pacific Decadal Oscillation (PDO) index was trending lower between 2000 and 2013, while the North Pacific index (NPI) and Antarctic Oscillation/Southern Annular Mode (AAO/SAM) index increased.

5 Results

We have identified several trends in the MISR cloud amount over the extratropical oceans in Fig. 1. Below, cloud fraction trends in each of the ocean basins shown in Fig. 1 are discussed in terms of the MCA decomposition for that basin. We proceed through each of the four basins and examine each MCA mode using shorthand NAX, NPX, SAX, and SPX to refer to MCA mode number X in the North Atlantic, North Pacific, South Atlantic, and South Pacific, respectively. Under each of the leftmost panels in Fig. 3, the percentage of covariance between the data sets and variance within each data set that can be explained by that mode are given as “(% covariance, % ERA-Interim variance, % MISR variance)”. Also, an estimate of the fraction of the basin-averaged total cloud trends, low cloud trend, and optically thick cloud trends (shown in Fig. 1) that can be explained by each MCA mode is given in Table 1. These values provide some additional information about the relative importance of each MCA mode.

5.1 North Atlantic

NA1 captures an extratropical trend of increased low-level cloud of moderate optical depth and reduced cloud of high optical depth, primarily at high altitudes. This mode is maximized along the storm track (roughly 40–55^∘ N) and is associated with increased temperature, pressure, and anti-cyclonic flow, as well as divergence at the surface, convergence aloft, and downward motion. Figure 2a–c show that trends in temperature, geopotential heights, and absolute vorticity in this region do not pass a 95 % confidence test; however, in the following discussion we see that similar MCA modes (modes with similar cloud, pressure, temperature, and vorticity patterns) exist in each of the other ocean basins studied where trends in the reanalysis variables are more robust. Figure 4a shows that there is no particularly strong connection between this mode and any climate index analyzed, but rather weak correlations with the East Atlantic teleconnection patterns and with the North Atlantic Oscillation. Li et al. (2014a) identify relationships between cloud occurrence and the NAO in CloudSat data; however, the correlations shown in Fig. 4 and the trends in the NAO and other East Atlantic teleconnection patterns are relatively weak for the time period analyzed. Nonetheless, NA1 describes 30 % of the basin average cloud occurrence trend, and the cloud occurrence pattern for this mode closely resembles the basin-averaged trend pattern in Fig. 1.

NA2 explains about 15 % of the covariance and 12 % of the variability in the ERA-Interim trends but only 4 % of the variance in the MISR trends. This mode is a smaller contributor to the overall (basin average) change in cloud occurrence (shown in Fig. 1), accounting for only 15 % of the basin-averaged trend compared with 30 % for NA1. The ERA-Interim profiles for NA2 show the largest trends in near surface q, T, and divergence, mid-troposphere vertical velocity, as well as the largest SST trend of any of the North Atlantic modes. While the changes in the cloud fraction joint histogram are somewhat noisy, the strongest response is in the low-level cloud. Interestingly, the spatial distribution to the left resembles recent trends in North Atlantic SST, with warming SSTs off the east coast of North America and cooling in the central North Atlantic (Fig. S1 in the Supplement). It has been documented that SSTs in the Gulf Stream influence cloud fraction (Minobe et al., 2008, 2010). Given the apparent connection to SST (and the fact that this mode does not have a strong correlation with any CPC index), we hypothesize that NA2 is a possible link between North Atlantic SST changes and low-level cloud fraction.

NA3 is primarily a subtropical mode (notice the location of red colors in left panel of Fig. 3) and is characterized by an increase in low-level cloud of moderate optical depth in the southern portion of the study region. This increase in low clouds appears to be due to an increase in specific humidity and increased convergence at low levels in the ERA-Interim data set, although as noted in Sect. 2, we caution that this increase in specific humidity may be a spurious feature in ERA-Interim data or may be overestimated and is not entirely corroborated by MERRA. A somewhat similar pattern of ERA-Interim and cloud fraction trends to NA3 is seen in the North Pacific in mode NP3, but NA3 features a small reduction in mid-level cloud around 4 km and larger reduction in high clouds above 7 km than NP3. This mode accounts for more of the total cloud occurrence trend in the basin than NA2 (18 %), accounts for 26 % of the low cloud trends, and indicates that the increase in low cloud seen in the Fig. 1 basin-averaged trends likely has separate contributions from the storm track region and the subtropical/equatorward part of the North Atlantic with different meteorological mechanisms.

5.2 North Pacific

NP1 is an east–west dipole, showing increased low-level cloud and decreased high-level cloud of median optical depth in the East Pacific, and the opposite to the west. The spatial structure of this mode resembles the spatial structure of the PDO, and Fig. 4 indicates that there is a noteworthy time-correlation between the two ($r=-\mathrm{0.42}$). It has been argued that the longer timescale PDO is a response to ENSO and other factors including variability in the Aleutian low (Schneider and Cornuelle, 2005; Newman et al., 2016). Newman et al. (2016) show that variability in the Aleutian low as captured by the NPI typically leads the PDO (see their Fig. 3b) by several months. Not surprisingly, we find an even stronger correlation between NP1 and the NPI (r=0.56) and to a lesser degree the PNA ($r=-\mathrm{0.48}$) (which is also known to be strongly correlated to the NPI). The panel on the right side of Fig. 3 indicates that NP1 cloud changes are associated with large changes in the thermodynamic variables near the surface, and a particularly large change in SST, with cooler near-surface temperatures and SSTs in the eastern third of the North Pacific associated with more low cloud, and conversely warmer SSTs being associated with less low cloud. We note that the sign of the MCA modes is arbitrary and, as shown, is opposite that defined by the PDO (PDO is positive with warmer SSTs in the eastern pacific and hence a negative correlation). NP1 also shows a weak increase in mid-tropospheric downwelling and low-level divergence with convergence aloft, which may explain the reduction in high cloud, and would further enhance stability in the lower troposphere. This mode, while it is the largest source of covariance between the cloud occurrence trends and the meteorological trends in this region, accounts for almost none of the cloud occurrence trend averaged over the full basin (2 % given in Table 1). This is because the spatial loading pattern for NP2 averages to near zero (meaning changes in the eastern and western part of the domain tend to cancel in the basin average). A useful feature of the MCA is that it can identify important patterns of cloud–climate interaction (strictly speaking covariability) that cannot be inferred from the information in Figs. 1 and 2 alone.

Most of trends for this region in Figs. 1 and 2 are captured by NP2, which is a very similar mode to NA1 (note NA1 not NP1) and describes 25 % of the total cloud occurrence trends in this region. It shows a reduction in high optically thick cloud and enhanced low cloud of moderate optical depth. The profiles of the ERA-Interim variables to the right in NP2 show increased temperature, pressure, and anticyclonic motion throughout the middle of the domain (along 45^∘ N), and are nearly identical to the NA1 pattern. Figure 2 shows that both 500 hPa temperature and pressure in this area have undergone robust positive trends. The time series derived from NP2 has a high correlation with the NPI (r=0.42), as does NP1, but NP2 has a much weaker correlation with the PDO than NP1. We examine the time series associated with these two MCA modes and the PNA, PDO, and NPI in Fig. 5a and b. These time series illustrate the close connection between both of the North Pacific MCA modes and the NPI, as well as the first mode's relation to the PDO and the second mode's relation to the PNA. However, given the difference in spatial patterns between NP1 and NP2, and the poor correlation of NP2 with the PDO, we interpret the pattern in NP1 as a response to longer timescale forcing related the PDO, whereas NP2 is associated with shorter timescale synoptic variability captured by the PNA mode and the NPI. In effect, while NP1 and NP2 both feature increased low clouds and decreased high clouds (albeit with some relatively subtle differences with regards to optical depth), they have very different spatial patterns which occur for different reasons.

Finally, NP3 shows a similar subtropical pattern to NA3, with an increase in low-level cloud of medium optical depth and an increase in low-level specific humidity mostly south of 40^∘ N. Like NA3 it contributes to the low cloud trend seen in the basin average (9 % in Table 1). This mode explains only a very small amount of the spatial variance in the cloud fraction data (2 %), but projects strongly onto the basin-averaged cloud occurrence trends (Table 1), which is again consistent with NA3.

Figure 5 Time series from the top two North and South Pacific MCA modes, shown along with the most correlated CPC indices identified in Fig. 3. The North Pacific index (NPI) is strongly connected to both North Pacific MCA modes (a, b), with correlations of 0.56 and 0.42 for modes 1 and 2, respectively. In the South Pacific, mode 1 (c) is strongly correlated with the Nino 4 index (−0.50), whereas mode 2 (d) is correlated with the Southern Annular Mode (0.32). In the plots above, the sign of each index that is negatively correlated with the MCA mode shown has been flipped and a 5-month boxcar filter has been applied to make the plots more readable. No filtering was applied to derive the correlation coefficients given (the correlations between all CPC indices and all MCA modes are given in Fig. 4a).

5.3 South Atlantic

The MCA results in the South Atlantic are the least tractable of any of the regions studied. SA2 appears to be a “correction” to SA1, with nearly identical large-scale spatial patterns, (but with SA2 having more fine-scale features), and nearly opposite changes in cloud fraction (except in the region between 1.5 and 2.5 km where SA1 shows little change and SA2 displays a decrease in the cloud amount). SA1 also explains much more of the variance in the MISR cloud fraction data than SA2 (21 % vs. 4 %), but explains very little of the spatially averaged trend. SA1 is associated with increased high pressure, temperature, and anti-cyclonic motion in the southeast part of the region (consistent with the overall ERA-Interim trends shown in Fig. 2g–i. The region of enhanced high pressure in SA1 corresponds with a region of high pressure in the SAM loading pattern, and indeed SA1 shows moderate correlation with the SAM (r=0.25) (Fig. 4a). SA3 explains only 7 % of the covariance and 1 % of the variance in the MISR data and is therefore not included in the analysis. In the NP and NA regions, we chose to include the third MCA modes because there was good agreement between the two basins and a coherent spatial pattern. There is no such agreement in the third MCA modes in the Southern Hemisphere regions, and the spatial patterns associated with the third modes appear noisy.

5.4 South Pacific

The changes captured by SP1 are related to the South Pacific Convergence Zone (SPCZ). The period studied was characterized by relatively neutral ENSO conditions, with more La Niña-like conditions later in the time series. The position of the SPCZ is influenced by ENSO, shifting to the south and west during La Niña conditions (Trenberth and Shea, 1987; Vincent, 1994), and the spatial pattern of cloud fraction changes shown in SP1 are consistent with the cloudy SPCZ shifting into the South Pacific study region. Figure 4a further corroborates the relation to ENSO, showing high time-correlation between this MCA mode and the Niño region 3, 3.4, and 4 anomalies ($r=-\mathrm{0.43}$ to −0.50). Figure 5c show the SP1 and the Niño 4 time series.

SP2 is most positive in the storm track region and the ECMWF profile trends shown on the right indicate enhanced geopotential heights, temperature, and anti-cyclonic flow, even displaying an increase in surface divergence with convergence aloft, similar to NA1, SA1, and NP2. These changes are associated with increased cloud of moderate optical depth at low levels and reduced mid- and high-level cloud of slightly larger optical depths. This mode explains 40 % of the basin-averaged cloud trends (Table 1). The SAM has undergone a positive trend in recent decades (Thompson and Wallace, 2000; Thompson et al., 2011, Fig. 3b), which is at least partly responsible for (or at least is correlated with) these changes. The sign convention is such that positive values of the SAM index correspond to an enhanced low over Antarctica and an intensified polar vortex, along with increased surface pressure and geopotential heights over much of the Southern Ocean. Figure 4a shows that the SAM has the highest time-correlation with SP2 (r=0.32), and the associated time series are shown in Fig. 5d. Arblaster and Meehl (2006) and Thompson et al. (2011) attribute most of the past trend in the SAM to anthropogenic ozone depletion, but increased CO₂ concentration (under an RCP 8.5 scenario, see Boucher et al., 2013) could cause this trend to continue even as Antarctic stratospheric ozone begins to recover (Thompson et al., 2011; Zheng et al., 2013). Hartmann and Ceppi 2014 relate changes in South Pacific reflected shortwave radiation observed by Clouds and Earths Radiant Energy System (CERES) Terra to a trend towards La Niña-like conditions and to trends in zonal mean winds and the SAM (which are very strongly correlated). However, they note that it is difficult to identify a robust relationship between the SAM and Southern Ocean cloud shortwave radiative effect trends due to the dominant influence of sea-ice changes and of ENSO over such a short time period. Ceppi and Hartmann (2015) note that while cloud amount, particularly mid- and high-level cloud, responds to changes in the annular modes, associated dynamically forced changes in cloud shortwave radiative effect may be of secondary importance to thermodynamically forced changes in the cloud phase. Regardless, the cloud amounts shown in SP2 and SA1 are likely driven by changes in the SAM, and continued data collection should help isolate cloud change responses to the SAM.

6 Conclusions

In closing, cloud data sets from EOS are now of sufficient quality and length to begin studying the response of cloud to synoptic variability on multi-year timescales. We have found a number of linkages between trends in synoptic meteorology and trends in cloud fraction via maximum covariance analysis. Notably, increased low cloud of moderate optical depth and reduced high cloud of higher optical depth are associated with increased temperature, anti-cyclonic motion, geopotential height, and subsidence in the extratropical storm track regions. We speculate that this could be linked to a strengthening of extratropical warm-core highs during the time period studied, although this would require additional analysis of daily data to verify. The maximum covariance analysis also revealed linkages between observed extratropical cloud fraction changes and known modes of climatic variability. Cloud changes (or responses) associated with ENSO, PDO, NPI, PNA, and SAM were found, as well as a possible link between cloud fraction and changing Atlantic SSTs (i.e. mode NA2). As climate records from MODIS and MISR (and we hope future generations of equally, if not more, capable instruments) continue, we will gain an increasing understanding of the response of clouds to synoptic variability on intra-decadal timescales. These observation-based relationships potentially offer exciting and new ways in which we can evaluate and improve climate models, as well as further our understanding of climate change.

Appendix A

Figure A1 Graphical depiction of how the maximum covariance analysis was formulated. Each of the boxes represents a two-dimensional (2-D) data matrix. The rows of each matrix correspond to latitude and longitude coordinates. Matrix “M” contains trends in the MISR cloud fraction joint histogram data. Matrix “E” contains trends in the ECMWF reanalysis data. The columns of “M” represent variability with respect to cloud optical depth and cloud-top height. The columns of “E” represent variability with respect to variable type and vertical level. Normalization (subtracting the mean and dividing by the standard deviation) is performed within each red box. C_ME is then a covariance matrix that captures spatial covariance between the trends in the cloud fraction data and the reanalysis data.

The MCA process is as follows: The MISR cloud fraction trend data were arranged into a two-dimensional matrix (M), with one dimension representing latitude and longitude, and the other representing cloud optical depth and cloud-top height trend for each MISR CTH-OD component. The ECMWF trend data were also arranged into a two-dimensional matrix (E). Latitude and Longitude were again represented by one dimension, whereas the trend variable type (temperature, geopotential, pressure, absolute vorticity etc.) at each height (z) comprise the second dimension. This is depicted in Fig. A1. Trends in each like-variable (e.g. Z, T, w, CF, and so on) were pre-normalized, meaning that trends for each variable had their mean value subtracted and were divided by their standard deviation. For instance, the mean of all temperature trends is subtracted from each temperature trend, and then the temperature trend data are divided by their standard deviation. A covariance matrix was then computed with respect to space, and a singular value decomposition was performed on this covariance matrix:

$$\begin{array}{}\text{(A1)}& \begin{array}{l}{\mathbf{U}}_{\mathbf{M}}\mathbf{\Sigma }{\mathbf{U}}_{\mathbf{E}}^{\mathbf{T}}={\displaystyle \frac{{\mathbf{ME}}^{\mathbf{T}}}{n}}.\end{array}\end{array}$$

Here, n is the number of latitude and longitude grid points included in the analysis. The right-hand side is a spatial covariance matrix computed between the two sets of normalized trends. U_M and U_E contain left and right eigenvectors produced by the singular value decomposition, where the columns of U_M and U_E represent modes of spatial covariance in MISR CTH-OD space and in the vertical profiles of the ECMWF reanalysis variables, respectively. The matrix Σ is diagonal, and includes the singular values associated with each mode. The spatial patterns in each trend data set associated with each column of U_M and U_E were retrieved by simply projecting U_M and U_E onto their respective data sets:

$$\begin{array}{}\text{(A2)}& {\mathbf{V}}_{\mathbf{M}}={\mathbf{U}}_{\mathbf{M}}^{\mathbf{T}}\mathbf{M}\text{and}{\mathbf{V}}_{\mathbf{E}}={\mathbf{U}}_{\mathbf{E}}^{\mathbf{T}}\mathbf{E}.\end{array}$$

Here, the columns of V_M and V_E are the spatial patterns in their respective trend data sets associated with each of the modes of covariance identified with the MCA. The columns of V_M and V_E were then standardized (in this case they are divided by their standard deviation, but the mean is not removed) and again projected onto the original dimensional versions of M and E (which we will call M^* and E^*). This yields dimensionalized versions of the MCA modes (${\mathbf{U}}_{\mathbf{M}}^{*}$ and ${\mathbf{U}}_{\mathbf{E}}^{*}$) and associated normalized spatial patterns (V_M and V_E), for example

$$\begin{array}{}\text{(A3)}& {\mathbf{U}}_{\mathbf{M}}^{*}={\mathbf{M}}^{*}{\mathbf{V}}_{\mathbf{M}}^{\mathbf{T}}\text{and}{\mathbf{U}}_{\mathbf{E}}^{*}={\mathbf{E}}^{*}{\mathbf{V}}_{\mathbf{E}}^{\mathbf{T}}.\end{array}$$

Each of the first two or three MCA modes in each ocean basin is displayed in Fig. 3. The leftmost panels show the spatial distribution associated with each mode derived using the ECMWF trends (V_E). These are normalized. The spatial distributions derived using the MISR data are omitted for space, but are necessarily quite similar, although not identical, to those shown. The middle panels show trends in cloud fraction joint histograms which are dimensionalized (${\mathbf{U}}_{\mathbf{M}}^{*}$), and the right panels show the associated dimensionalized trends in the profiles of the various ECMWF reanalysis variables (${\mathbf{U}}_{\mathbf{E}}^{*}$). Recall, that the true trend in the original data set associated with each mode can be recovered by projection of the contents of the middle or right panel (${\mathbf{U}}_{\mathbf{M}}^{*}$ or ${\mathbf{U}}_{\mathbf{E}}^{*}$) onto the spatial distribution in the left panel (V_E). Changing the sign of all three panels would then yield the same result, and we have chosen the sign of each mode such that they most resemble the patterns shown in Fig. 1.

In ascertaining the significance of each MCA mode, it is important to examine the covariance between the two data sets and the variance within each data set explained by the MCA mode. These values are printed in Fig. 3 underneath each spatial loading pattern in the following format: “(% covariance, % ERA-Interim variance, % MISR variance)”. These measures give an idea of a mode's relative importance in the context of the original data. The percent covariance explained is a measure the importance of a mode relative to other modes, whereas the percent ERA or MISR variance explained measures the importance of a mode for explaining the trends in each particular data set. For instance, the first and second modes in the North Pacific are quite robust because they explain a large amount of the total covariance between the two data sets and explain a substantial amount of the variance in each data set. Conversely, the third North Pacific mode is less important, because while it explains a non-negligible fraction of the total covariance, it explains only 2 % of the variance in the MISR data and thus does not project strongly onto the MISR trends. While there is no universally agreed upon method for testing the significance of MCA results, a useful metric is the “normalized root mean squared covariance” (RMSC):

$$\begin{array}{}\text{(A4)}& \mathrm{RMSC}={\displaystyle \frac{{\left|\left|{\mathbf{ME}}^{\mathbf{T}}\right|\right|}_{\mathrm{F}}}{\sqrt{\mathrm{tr}\left({\mathbf{MM}}^{\mathbf{T}}\right)\mathrm{tr}\left({\mathbf{EE}}^{\mathbf{T}}\right)}}}.\end{array}$$

Here, $\left|\right||{|}_{\mathrm{F}}$ is the Frobenius norm and tr() is the trace operator. In general, larger values of this metric imply a robust result whereas smaller values imply that the MCA modes do not capture a significant portion of the variance in the two data sets. For artificially generated data with similar dimensions to the data used here, poorly correlated fields yield RMSC ≈0.09 whereas well correlated fields yield RMSC ≈0.26, (this range will vary depending on the size of the data set and the number of independent samples). The RMSCs computed for each of the study regions are as follows: North Atlantic, 0.17; North Pacific, 0.17; South Atlantic, 0.16; and South Pacific, 0.17. These values are reassuringly high considering the large variability in the monthly MISR data. Finally, each of the MCA modes, which are shown in Fig. 3, indicate the percent covariance explained, and the percent variance explained in each of the two data sets by that mode. We discuss these modes in Sect. 5.

In Table 1, we provide an estimate of the fraction of the cloud trends in Fig. 1 that can be explained by the cloud–meteorology relationship represented by each MCA mode. This is computed by projecting the MCA mode's cloud occurrence pattern (${\mathbf{U}}_{\mathbf{M}}^{*}$) onto the spatial loading pattern derived from the meteorological data (V_E) to get an estimate of the cloud fraction trend represented by that mode (Δ_MCA). The fraction of the Fig. 1 trend explained is then calculated as follows:

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\sum \left|{\mathrm{\Delta }}_{\mathrm{MCA}}-{\mathrm{\Delta }}_{\mathrm{obs}}\right|}{\sum \left|{\mathrm{\Delta }}_{\mathrm{obs}}\right|}},\end{array}$$

where the summation is over each of the bins in the cloud occurrence joint histograms. By applying the summations in Eq. (5) to subsets of the cloud occurrence histogram cloud categories we also obtain estimates of the fraction of the low (CTH < 2.5 km) and thick (optical depth > 23) cloud trend explained. Additionally, the combination of the MISR cloud fraction joint histogram trends and associated spatial distribution (U_M and V_M) can be projected back onto the original cloud fraction data set to yield the time series associated with each of the MCA modes. More precisely, this can be done by taking the outer product of a column of V_M and a (transposed) column of U_M, which results in a matrix of the same size as M. Each month of the original cloud fraction data can be restructured into a matrix of the same size as M, and then the sum of the element-wise product of the outer product of the column of U_M and V_M and each month of MISR cloud fraction data can be taken to retrieve a time series. This time series is then standardized by subtracting its mean and dividing by its standard deviation. Each of the resulting time series based on the MCA modes will necessarily indicate trends because the MCA was performed on the trends computed using the original data sets.