Large-scale conditions over subtropical marine stratocumulus areas are extracted from global climate models (GCMs) participating in phase 3 of the Coupled Model Intercomparison Project (CMIP3) and used to drive an atmospheric mixed-layer model (MLM) for current and future climate scenarios. Cloud fraction is computed as the fraction of days where GCM forcings produce a cloudy equilibrium MLM state. This model is a good predictor of cloud fraction and its temporal variations on time scales longer than 1 week but overpredicts liquid water path and entrainment.
GCM cloud fraction compares poorly with observations of mean state, variability, and correlation with estimated inversion strength (EIS). MLM cloud fraction driven by these same GCMs, however, agrees well with observations, suggesting that poor GCM low cloud fraction is due to deficiencies in cloud parameterizations rather than large-scale conditions. However, replacing the various GCM cloud parameterizations with a single physics package (the MLM) does not reduce intermodel spread in low-cloud feedback because the MLM is more sensitive than the GCMs to existent intermodel variations in large-scale forcing. This suggests that improving GCM low cloud physics will not by itself reduce intermodel spread in predicted stratocumulus cloud feedback.
Differences in EIS and EIS change between GCMs are found to be a good predictor of current-climate MLM cloud amount and future cloud change. CMIP3 GCMs predict a robust increase of 0.5–1 K in EIS over the next century, resulting in a 2.3%–4.5% increase in MLM cloudiness. If EIS increases are real, subtropical stratocumulus may damp global warming in a way not captured by the GCMs studied.
Low clouds reflect solar radiation, strongly cooling the planet. Even small changes in the prevalence of these clouds would induce an energy imbalance of magnitude similar to that of CO2 doubling (Randall et al. 1984; Slingo 1990; Hanson 1991), making low cloud feedback one of the largest sources of climate change uncertainty (Bony and Dufresne 2005; Solomon et al. 2007). Low clouds are particularly difficult to simulate in global climate models (GCMs) because they are typically thinner in the vertical than a single GCM grid cell and result from an intricate balance between complicated microscale processes.
In light of the importance of low cloud feedback, much effort has gone into devising clever ways to circumvent issues with GCM low cloud predictions. For example, Medeiros et al. (2008) grouped aspects of cloud change from two GCMs run in aquaplanet and earthlike configurations into model-dependent and model-independent categories. Thermodynamic effects were found to be more important than circulation changes in explaining cloud changes (a result corroborated by Bony and Dufresne 2005). Clement et al. (2009) ranked GCMs from the phase 3 of the World Climate Research Program’s Coupled Model Intercomparison Project (CMIP3) according to their ability to reproduce current-climate relationships between low clouds and their environment. Their highest-ranked model [the Met Office (UKMO) Hadley Centre Global Environmental Model version 1 (HadGEM1)] predicts a reduction in low clouds. X. Qu et al. (2012, manuscript in preparation) explore the connection between current-climate variability and climate change projections for low clouds in the CMIP3 models. They find that models with stronger sensitivity to estimated inversion strength (EIS) (Wood and Bretherton 2006, hereafter WB06) in the current climate have more positive low-cloud fraction change (negative cloud feedback) under future warming but that model-predicted ΔEIS is a poor predictor of the magnitude of future change. Eitzen et al. (2011) break satellite-observed low cloud variability into terms correlated with reanalysis-based subsidence ws and EIS. Arguing that EIS and ws should be relatively unaffected by climate change, they use the component of cloud response to SST uncorrelated with EIS or ws to diagnose a positive low cloud feedback. Satellite and surface-observer networks have also been used to diagnose trends, though results vary depending on dataset (Norris 1999, 2005) and concerns about data quality exist (Evan et al. 2007).
Another tactic for understanding low cloud change is to examine results from regional models driven by large-scale conditions for current and future climate. Simulating a limited portion of the globe allows for higher resolution and more specialized low-cloud parameterization but requires specification of the large-scale flow into the region of interest. Most studies of this type have used idealized boundary conditions following the two-box tropical model popularized by Pierrehumbert (1995) but with cold-pool to warm-pool coupling ignored. Zhang and Bretherton (2008) used this approach with a single-column model (SCM) to explore the reasons for negative low cloud feedback in the Community Atmosphere Model version 3, finding low cloud feedback to result from unphysical interaction between parameterizations of shallow convection, deep convection, and boundary layer (BL) turbulence. Caldwell and Bretherton (2009, hereafter CB09) explored low cloud feedback using an atmospheric mixed-layer model (MLM). They identified a new “subsidence–lapse rate feedback,” which increases cloud water due to a combination of decreased subsidence and increased lapse rate, causing the BL to deepen while simultaneously reducing entrainment drying. Xu et al. (2010) found negative low-cloud feedback in cloud resolving model (CRM) simulations when SST is increased by 2 K and forcings were idealized following Zhang and Bretherton (2008). Stratus increase occurred because the BL deepened while remaining well mixed. Shallow cumulus increase, which was especially strong, was due to a rise in turbulent moisture transport. Blossey et al. (2009) used output from a GCM whose cloud physics were replaced by a 2D CRM to drive higher-resolution CRM simulations. Whereas Wyant et al. (2009) found cloudiness in the parent GCM to increase due to stronger cloud-top cooling in a warmer and moister BL, the higher-resolution simulations of Blossey et al. had more muted cloud increases, which were mainly due to increased inversion strength. Lauer et al. (2010) ran a regional model forced by reanalysis with and without shifting by GCM-predicted climatological forcing change (a “pseudo-global warming” study). In contrast to other modeling results, this group found low cloud reductions in a warmer world due to decreased BL depth, a result they attributed to increased free-tropospheric opacity in a warmer climate. Rieck et al. (2012) also found positive shallow cumulus cloud feedback in large-eddy simulations (LES) with fixed background relative humidity due to BL deepening and drying in response to increased surface fluxes.
In summary, GCMs generally predict future decreases in low clouds, though intermodel spread is large and simulating these clouds continues to be a challenge. Low clouds seem to be more sensitive to thermodynamic rather than circulation changes, suggesting that improving GCM cloud physics parameterizations is the best way to reduce intermodel spread. Local-area modeling studies tend to show increased low clouds with warming, though the strength and mechanisms for change seem to vary between models. It is still unclear how much of the differences between regional modeling studies come from model formulation details and how much comes from boundary-condition assumptions.
In this study, we follow the regional modeling approach by driving a simple local-area model with forcings from all available CMIP3 models. Our single-model, multiforcing approach nicely complements ongoing work by the Cloud Feedbacks Model Intercomparison Project (CFMIP)–Global Energy and Water Cycle Experiment (GEWEX) Cloud System Study (GCSS) Intercomparison of Large Eddy Models and Single Column Models (CGILS) project (Zhang et al. 2010, 2012, manuscript submitted to J. Adv. Model. Earth Syst.; Blossey et al. 2012, manuscript submitted to J. Adv. Model. Earth Syst.), which aims to drive many regional models with a single set of forcing conditions. While CGILS will provide a sense of the uncertainty provided by local physics assumptions, our approach will improve understanding of model sensitivity to large-scale forcing specifications as well as clarifying how large-scale forcings and their changes should be parameterized. Both our study and CGILS aim to quantify the relative contributions of large-scale conditions and local physics to low-cloud feedback uncertainty by ignoring uncertainty from one contributor or the other. Because local and large-scale factors may interact nonlinearly,1 both approaches are needed to reliably partition the source of intermodel spread. Additional motivation for our approach is that the simplicity of our model allows us to give physical explanations for our findings and using an ensemble of forcing conditions provides a sense of forcing uncertainty that has been absent in previous single-model/single-forcing studies.
In the next section, we explain the details of our methodology. In section 3, we test the skill of our framework at predicting present-day stratocumulus (Sc). Section 4 focuses on CMIP3 model predictions of large-scale forcing changes. The effect of these changes on MLM clouds is discussed in section 5. Section 6 provides a closer look at why cloudiness changes in the future and conclusions follow in section 7.
2. Data and methods
We use daily data from 10 models in the CMIP3 archive (Meehl et al. 2007) to investigate changes to the large-scale environment in which subtropical Sc are embedded. We use a mixed-layer model forced by this data to translate forcing changes into impacts on cloudiness. Our modeling approach is diagrammed in Fig. 1. Models were chosen based on availability of daily forcings needed to drive the MLM. We use years 1981–2000 from coupled twentieth-century runs as “current conditions” and years 2081–2100 from scenario A1B as “future conditions.” Table 1 lists the models used. More complete documentation is available online (at http://www-pcmdi.llnl.gov/ipcc/model_documentation/ipcc_model_documentation.php). For validation purposes we also force our MLM with data from the European Centre for Medium-Range Weather Forecasts Interim Re-Analysis (ERA-Interim) (Dee et al. 2011). We compare the resulting cloud fraction against International Satellite Cloud Climatology Project (ISCCP) D1 data daytime mean cloud fraction with cloud-top pressure lower than 560 mb2 and cloud optical thickness greater than 3.55 (Rossow and Schiffer 1999). Because our model is only appropriate at the core of Sc areas, our study focuses on the five major subtropical marine Sc regions identified in Klein and Hartmann (1993, hereafter KH93). These regions are listed in Table 2. For each GCM, all grid cells lying within each study region are used independently to run the MLM to equilibrium for each day of the experiment.
We choose a MLM as our local-area model because it is straightforward to interpret and efficient enough to use for multiple multidecade simulations, yet it is also an accurate model of the Sc-topped BL as long as turbulence is strong enough to keep the BL well mixed. Because these conditions are definitely not met for shallow cumulus, our study does not explore the cumulus component of low cloud feedback, which Medeiros et al. (2008), Wyant et al. (2009), and Webb et al. (2013) show is also important. Lilly (1968) developed the MLM approach by observing that moist- adiabatically and dry-adiabatically conserved variables tend to be spatially uniform in the Sc-topped BL, allowing him to collapse the 3D equations for Sc-topped BL evolution into scalar budget equations for BL energy, mass, and moisture. Detailed specification of these budget equations is given in Zhang et al. (2005). Entrainment (we), the incorporation of warm, dry free-tropospheric air into the BL by cloud-top turbulence, is parameterized in our model following Lewellen and Lewellen (1998) with large-eddy we efficiency η = 0.25. This parameterization was chosen because it reproduces observations well with our default η value (Stevens et al. 2002; Caldwell et al. 2005) and can be made to mimic other parameterizations by varying η (Zhang et al. 2005). Sensitivity to η is explored in appendix B. Radiation is included following the simple empirical formulation of Stevens et al. (2005b) with current-climate maximum cloud-top cooling of 40 W m−2 chosen to mimic diurnally averaged conditions as in Zhang et al. (2009, hereafter Z09). We reduce this number by 4 W m−2 for future climate conditions based on offline radiative transfer calculations using global versions of the Iacono et al. (2008) schemes with CO2 doubled to match the year 2000–2100 CO2 change in the A1B scenario. Divergence D is assumed constant with height so ws = Dz where z is altitude in meters. Following Zhang et al. (2005), drizzle is parameterized as proportional to the cube of cloud depth with subcloud evaporation ignored and droplet concentration held fixed at 100 cm−3. Horizontal advection of BL temperature (moisture) is parameterized using the (saturation mixing ratio of) SST gradient.
Free-tropospheric, surface, and horizontal-advective forcings from GCMs are needed by the MLM. Most of these are available directly from the CMIP3 archive but a few must be derived from other data as described in appendix A. Where possible, forcings computed using these derivations are tested against output directly from reanalysis and found to add little uncertainty. A critical issue is that free-tropospheric properties are required as boundary conditions for the MLM, but varying BL depth makes it impossible to choose a particular model layer as always being the bottom of the free troposphere. This problem is illustrated in Fig. 2, which shows ERA-Interim relative humidity (RH) over the California Sc region between December 2006 and March 2007. We create free-tropospheric forcings by using 850-hPa values for spatial locations and times where the bottom of the free troposphere is diagnosed as lower than 850 mb. When 850 hPa is in the BL, we reconstruct synthetic free-tropospheric 850-hPa values by extrapolating free-tropospheric levels downward. Levels diagnosed as the bottom of the free troposphere are included as a white line in Fig. 2. The details of this process are included in appendix A. Model sensitivity to this correction is shown in appendix B.
A limitation to conventional mixed-layer models is that horizontal homogeneity forces them to have cloud fractions of zero or one. Z09 extends the MLM theory to predict cloud fraction by running the MLM to equilibrium for a series of observed conditions and taking cloud fraction to be the fraction of time that the MLM predicts a cloudy solution. Stratocumulus-free conditions are said to occur when the MLM predicts a clear solution, when the BL grows more than 2 km deep (at which point the BL is assumed decoupled into cumulus), when equilibrium is not reached in 200 days (which almost never happens), or when D < 0 m−1. This last criterion is just to avoid needless simulations—equilibrium BL depth will almost certainly be above 2 km when mean motion is upward. The relative importance of these criteria on Sc incidence is explored in section 6. We run our model to equilibrium for every day of the simulation period using diurnally averaged large-scale forcings. As noted in Z09, longer averaging periods are problematic because divergence has a lot of power at short time scales and cloud fraction decreases nonlinearly as divergence becomes washed out by averaging. Each simulation is initialized from cloudy conditions based loosely on data from the Second Dynamics and Chemistry of the Marine Stratocumulus field study (DYCOMS-II) (Stevens et al. 2007). As shown in Fig. 5 of Z09, initial conditions do make a difference to the equilibrium solution (because both cloudy and clear equilibria are possible from the same MLM forcing). We choose to initialize from cloudy conditions because upstream conditions in Sc regions are typically from areas of greater cloudiness. MLM forecast skill is presented in section 3a. In general, the Z09 method works well in the regions studied here since conditions change slowly in the upstream direction; in other regions, such as near the coast or equator, cloud fraction is overpredicted.
The MLM cloud fraction approach was shown in Z09 to successfully capture many observed features of low clouds. Since our model differs from Z09 by including drizzle, not assuming 850 hPa is always in the free troposphere, and replacing surface divergence with a value chosen to approximate 850-hPa ws, we should confirm that it still performs well. Since the current study also extends upon Z09 by using CMIP3 forcings instead of reanalysis, we should also verify that MLM runs driven by CMIP3 forcings lead to a good MLM simulation of current-climate Sc properties.
a. Reanalysis-forced runs
In the current climate, variations of Sc cloud fraction are very well correlated with variations in EIS (WB06), so we consider reproducing this relation in our MLM to be of critical importance. Figure 3 shows the relation between EIS and cloud fraction based on ERA-Interim data from 2006 (a recent year of weak ENSO amplitude). As in Z09, results agree very well with ISCCP observations. The correlation between seasonal- and area-mean cloud fraction and EIS in the five dominant Sc regions is 0.91, which is very similar to the 2006 ISCCP/ERA-Interim value of 0.93 and the WB06 value of 0.92 for all KH93 regions. The slope of the MLM regression line is 8.7% K−1, somewhat higher than the ISCCP observed value of 6.9% K−1 for year 2006 and the WB06 value of 6.0% K−1 from climatological surface-observer data. Since our model overpredicts cloud sensitivity to EIS changes, its predicted cloud changes may be amplified. This potential for bias should be kept in mind but will not affect our conclusions since we seek to understand how and why clouds will change rather than to make quantitative predictions. The annual-mean cloud fraction from the MLM in the five regions is 5.8% higher than found for ISCCP, though discrepancy between daytime-average observations and diurnally averaged model output may explain some of this difference. MLM cloud fraction predicted here is also larger than found in Z09; this is due to the addition of drizzle, which increases cloud fraction by decreasing the frequency of cases where BL depth is above 2 km, the height beyond which the BL is assumed decoupled. The fact that good model skill is found even with altered parameterizations and a different forcing dataset fosters confidence in the Z09 methodology.
Intuitively, the Z09 approach works because the MLM predicts cloud whenever the driving large-scale conditions are conducive to cloud and hence cloud is likely observed to be present. This means that the temporal evolution of MLM cloudiness should track that of observed cloudiness. This is explored in Fig. 4, which shows the correlation between MLM and ISCCP cloud fraction after averaging to different time scales. In all regions, MLM forecast skill increases for longer time scales. Forecast skill also shows regional differences. In the Peru region, MLM-predicted cloud fraction is highly correlated with observations, even at a daily time scale. In the Canary region, MLM and ISCCP cloud are poorly correlated at all but the longest time scales. Except for the Canary region, the square of the correlation coefficient is greater than 0.5 for all time scales longer than a week. Prevailing forecast scores such as the Brier skill score (Brier 1950; Jakob et al. 2004) show similar results (not shown).
There are several possible explanations for poor MLM skill at short time scales. Our study associates each day’s forcing with its equilibrium cloud state, while real clouds may take a day or more to respond to a forcing (Schubert et al. 1979). A mismatch of a couple of days is not a problem on the climatological scales we are focused on but ruins correlation at short scales. Additionally, while small-scale forcing variations in space and time can radically influence instantaneous cloud amounts [e.g., in the pockets of open cells explored by Stevens et al. (2005a) and Comstock et al. (2005)] these fluctuations are not captured by our daily, gridcell-mean forcings. When averaged over several days, however, the influence of variability on short time and space scales is diminished, improving MLM skill. In addition, random instantaneous errors in ERA-Interim or ISCCP observations (Marchand et al. 2010) may contribute to poor correlation at shorter time scales that disappear after averaging. The above explanations suggest that bad short-term forecasts do not imply poor climate projections. Additionally, our model is extremely simplified, so perfect agreement is not expected. In particular, we predict clear sky whenever our crude decoupling threshold is reached, while in reality decoupled Sc decays slowly into patches of cumulus. Using equilibria, omitting downstream development, and ignoring the diurnal cycle could also cause problems. In light of these simplifications, the level of agreement with observations is quite impressive. In summary, current-climate MLM results are sufficiently accurate that MLM climate change predictions should be taken seriously.
b. Twentieth-century CMIP3-forced runs
One test of the value of driving a MLM with GCM forcing is whether current-climate statistics improve relative to those from the parent GCM. Since geographic and seasonal variability in cloud fraction is strongly correlated to EIS variations in the current climate, the relationship between these variables provides a compact test of model quality. The left panel of Fig. 5 shows that GCMs generally fail to reproduce observed seasonal and geographical response to EIS change (as represented by the WB06 line), instead predicting ~60% cloud fraction most of the time. This is apparently not due to bad large-scale conditions since when the MLM is driven by these same large-scale conditions (Fig. 5b) observed patterns of variability re-emerge. This is strong evidence that poor Sc prediction in GCMs is largely due to grid-scale parameterizations (physics) rather than large-scale forcing (dynamics). The large-scale information needed for predicting cloud fraction is present but not adequately used by the physical parameterizations of the CMIP3 models. GCM insensitivity to EIS is important since EIS is shown in section 5 to play a critical role in MLM response to warming.
One caveat to the above analysis is that total cloud instead of low cloud fraction is used for the GCMs because low cloud fraction is only calculable for 4 of the 10 models that have the needed MLM forcings. For these four models, we computed GCM low cloud fraction assuming random overlap of daily resolution 3D cloud data below 700 mb. Mean differences and correlations between total-cloud and random-overlap low-cloud fraction derived in this way are generally quite good (Table 3). Because higher clouds are not generally observed in the core Sc regions, this result is unsurprising. We note, however, that this relationship may not hold for all models; in particular, Broccoli and Klein (2010) found low-cloud fraction in the GFDL CM2.1 to respond very differently than total cloud to current-climate meteorological variations. This may be due to their use of a larger study region (15°–25°N, 115°–145°W) that may contain more mid- and upper-level clouds than our corresponding northeast Pacific (California) region. In any case, our approximation makes physical sense and seems to hold in the regions we have defined for the models we can check, justifying our cautious use of GCM total-cloud fraction for the remainder of the study.
The MLM-predicted mean BL structure for cloudy periods is summarized in Table 2. Comparison against satellite-derived climatological values (included in parentheses in Table 2 where available) suggests in-cloud liquid water path (LWP) is substantially too high. LWP overprediction is expected because of the use of adiabatic assumptions and neglect of stratification. The equilibrium assumption likely adds to overprediction by precluding lower-LWP clouds transitioning between clear and fully thickened equilibrium states. Entrainment is also seen to be overpredicted, though comparison against observations is more tenuous because satellite-derived entrainment is only available for the season of maximum Sc and is sampled across all conditions instead of just when Sc are present. Modeled BL depth appears to be reasonable, though again the satellite data includes times without Sc, which might bias the observations high. Ship-based data suggest BL depths substantially lower than found for the MLM (Bretherton et al. 2004; Kollias et al. 2004; Albrecht et al. 1988; Stevens et al. 2003) but tend to focus on the season of maximum Sc extent, when BL depths are lower than average. In all, model results are fairly reasonable given the simplicity of the MLM framework.
4. Large-scale forcing changes
The forcings needed by the MLM form a comprehensive list of the large-scale conditions that significantly impact low clouds; studying them directly can provide insight into low cloud response to climate change without requiring our modeling assumptions. Table 4 shows the twentieth-century changes in MLM forcings from the CMIP3 models. Variables where at least 8 of 10 models have changes of the same sign are indicated by boldface type. Several of these changes are expected from previous research. Warming of the sea surface and at 850 hPa is a foregone conclusion for a global warming scenario. Increased dθ/dz is the expected result of lapse rate feedback (Hansen et al. 1984). Divergence decreases as expected from Knutson and Manabe (1995), Held and Soden (2006), and Vecchi et al. (2006). Surface wind speed decreases slightly, yet cold advection strengthens because of increased SST gradients. Explicit use of the Clausius–Clapeyron relation in constructing υ · ∇q ensures increased dry advection in a warmer climate. Free-tropospheric RH changes are weak and inconsistent as expected from Solomon et al. (2007) and Soden et al. (2002); free-tropospheric q increases with T as predicted by the Clausius–Clapeyron relation.
An unexpected result is that EIS increases in all regions by >15% relative to current-climate values. Expressed in terms of degrees EIS change per degree global-averaged surface warming, this increase is between 0.3 and 0.5 K K−1, which is consistent with the results for the 80%–100% lower-tropospheric stability bin in the CMIP3 slab ocean runs analyzed in Webb et al. (2013). When developed, EIS was expected to be relatively invariant to climate change (WB06). As a result, Eitzen et al. (2011) use this assumption to separate a climate change signal from internal variability in satellite data. EIS change is also assumed to be weak in CGILS specification. Our results suggest that constant EIS is a bad assumption. EIS change is shown in section 5 to have a big effect on MLM low clouds so understanding why EIS is changing is important. Since free-tropospheric temperature in the subtropics is influenced by convection in the warm pool and surface temperature over the ocean is strongly coupled to local SST, EIS increase may be a manifestation of increased warm pool/cold pool SST gradient (CB09). Increased EIS could also be the result of increased near-coast subsidence warming and warmer free-tropospheric advection as land–sea temperature gradients are predicted to increase in a warming world (Sutton et al. 2007). These hypotheses will be explored in a future paper.
To illustrate the importance of correcting for incursions of the BL beyond 850 mb, RH taken directly at 850 hPa is also included in Table 4. Since this table also shows that CMIP3 BL depth (diagnosed using the algorithm described in appendix A) increases in a warmer climate, one might expect incursions of moist BL air at 850 mb to increase, raising RH. Surprisingly, Table 4 instead shows drying at 850 mb. This seems to be due to decreased RH during times when 850 hPa is in the BL (which happens >50% of the time in CMIP3 model simulations of current climate).
Nonlinearity of clouds means that changes in variability could also affect GCM low-cloud change. For example, Sc is reduced under very weak or very strong subsidence, so narrowing the width of the divergence probability distribution function (PDF) without changing its mean value may increase cloudiness (Z09). Changes in the temporal standard deviation σ of the MLM forcings averaged over all cells in each region are presented in Table 5. The seasonal cycle is removed before σ is computed by subtracting a three-point running-mean filtered daily resolution climatological seasonal cycle. Except for temperature advection (which strengthens in the mean but weakens in variability) and BL depth (which weakens in variability, probably because archive grid spacing increases with height), changes in mean and variability tend to be of similar sign though the relative strength of mean versus σ changes differ between variables. For example, divergence fluctuations decrease significantly more than mean changes in most regions. EIS variability, on the other hand, has little change even though EIS increases substantially in the mean. Variability changes for many variables are consistent across models and across regions, suggesting that physical reasons underlie these model changes. Variability changes do not seem to have a systematic effect on mean MLM cloud change, however, since running our MLM using current-climate forcings translated to have mean values that match the future scenario (i.e., pseudo-global warming runs) produced results very similar to those described below (not shown).
A goal of this work is to identify not only the extent to which the various large-scale forcings are consistently predicted by the CMIP3 ensemble but also the degree to which forcing differences influence cloud change predictions. This task is made more complicated by the fact that large-scale forcings are strongly interconnected. These interconnections are highlighted in Fig. 6, which shows monthly average3 temporal correlations between forcing variables (after removing climatological annual cycles). For each set of forcings, correlations are given separately for each model and for each region. Pairs of variables that show consistent colors across models and regions indicate relationships that are captured consistently across the CMIP3 archive. Large divergence is seen to be associated with a strong inversion, a dry free troposphere, and increased surface winds with resulting increased cold, dry advection. These relationships were also found in the Klein et al. (1995) analysis of interannual variability of the summertime subtropical high in the northeast Pacific. Temperature and moisture advection are perfectly correlated owing to their construction from the same SST gradients and surface winds. This web of correlation should be kept in mind when analyzing model response to forcings as cloud change may be correlated to a forcing only through association with other variables.
A potential problem with the regional model methodology is that a driving model may imprint its own cloudiness into its forcings. A possible way this could happen is through less low cloud leading to increased insolation, which would cause increased SST, decreased EIS, and hence less MLM cloudiness. Fortunately, however, Fig. 6 generally shows little correlation between GCM cloud fraction and forcing variables, suggesting that temperature imprinting is not important. One exception is a modest negative correlation between GCM cloud fraction and divergence, perhaps because CMIP3 GCMs with a deeper BL are more likely to support clouds. This is the opposite of the relationship found in the MLM, where increased subsidence decreases the likelihood of MLM decoupling and therefore increases cloudiness. Further evidence that imprinting is not a problem is given in section 5, where MLM and GCM cloud fraction are shown to be uncorrelated.
5. Climate change
The forcing changes documented in Tables 4 and 5 have complex effects on stratocumulus. Weakened divergence allows for a higher BL top for a given we. This could increase clouds because a higher and hence colder BL top is more likely to be saturated. It could also decrease cloudiness because deeper BLs are more likely to be decoupled. Increased EIS weakens entrainment, increasing cloudiness by increasing BL relative humidity and by reducing the likelihood of decoupling. Free-tropospheric moisture increases result in decreased evaporation of entrained air. This increases cloudiness by reducing evaporation and, potentially, (though not in our MLM) by reducing evaporative generation of turbulent entrainment. Free-tropospheric moisture can reduce turbulence by blocking cloud-top longwave cooling, but we do not include this effect in our model because offline radiative transfer calculations suggest that this radiative damping is canceled by increased radiative cooling from concomitant BL warming and moistening. Increased cold and dry advection have opposing effects on condensation level and hence cloudiness. Decreased surface winds suggest decreased surface fluxes, which again have opposing effects on cloud. In this section, MLM runs are used to predict cloud changes.
In Fig. 7 the MLM climate change signal in cloud fraction for each CMIP3 GCM is plotted against cloud change from its host GCM. As noted earlier, the direct impact of carbon dioxide on infrared opacity is included by reducing maximum cloud-top radiative cooling from 40 to 36 W m−2. Smaller points show what results would have been if this effect was neglected; the direct radiative impact is negligibly small.
Data for each model are normalized by that model’s global-average temperature change to avoid convolving cloud sensitivity with magnitude of forcing change. This normalization reduces spread in GCM but not MLM cloud fraction (not shown). One possible explanation for this is that GCM clouds are sensitive to the magnitude of SST change while MLM clouds are more sensitive to the pattern of change. Alternatively, differences in global temperature change in the GCMs could be a response to low cloud changes. Our MLM framework precludes this feedback because Sc are not allowed to affect the large scale. MLM and GCM cloud changes are poorly correlated (r2 = 0.14), reinforcing our conclusion that GCMs are not imprinting their cloudiness onto the MLM.
Surprisingly, filtering GCM results through a MLM does not reduce intermodel spread after normalizing. The standard deviation of MLM change across models is 1.3% K−1 compared to 1.1% K−1 for the GCMs. As seen in section 3b, MLM cloud fraction is much more sensitive to EIS; diversity in this variable has little impact on GCM cloud but effectively disperses the MLM results. This suggests that improving cloud physics parameterizations is not a sufficient condition for reducing intermodel spread in low cloud change. Important variations in large-scale forcings exist; CMIP3 models just are not sensitive to them (a fact also pointed out in Webb et al. 2013).
While CMIP3 GCM future cloudiness decreases by 0.8% K−1 on average, cloudiness change averaged over the MLMs actually increases by 0.9% K−1. All data points in Fig. 7 lie at or above the 1:1 line, indicating that all GCM physics parameterizations produce more negative cloud changes than predicted by MLM physics. The 1.7% difference between MLM- and GCM-predicted low cloud change can be put in context by noting that a 1% increase in low cloud induces a 1 W m−2 local decrease in absorbed net top-of-atmosphere (TOA) radiation (KH93), and ~10% of the planet is covered by subtropical oceanic Sc,4 leading to a global radiative impact of about −0.17 W m−2 K−1. Since the global net cloud feedback is thought to be around 0.5 W m−2 K−1 (Soden and Held 2006), the cloud changes found here are very noteworthy.
The source of differing MLM and GCM cloud response is explored in Fig. 8, which shows the correlation between forcing and MLM cloud fraction for each Sc region at three different time scales. Figure 8a shows the correlation across models of cloud change versus forcing change between the late twentieth century and the twenty-first century. Large values in this panel indicate that the relative magnitude of MLM cloud change can be skillfully predicted from its relative change in that forcing. Figure 8b shows correlation across models of current-climate mean forcing versus mean cloudiness. Because our study is based on 10 CMIP3 models, correlations in Figs. 8a,b are based on 10 data points each. Figure 8c gives temporal correlations between cloud fraction and large-scale forcings within each model using monthly average data for years 1981–2000 with the seasonal cycle removed (totaling 19 × 12 time samples per calculation). Cells in this panel are subdivided by model in a similar manner as in Fig. 6. A striking feature of Fig. 8 is that EIS is a great predictor not only of current-climate temporal variability within a model but also of the mean current-climate MLM cloud fraction a given model will have and (to a lesser extent) how that model’s low cloud will change in the future. Therefore, it is reasonable to explain MLM predictions of future cloud increases as a consequence of the large increases in estimated inversion strength found throughout the CMIP3 ensemble. The physical basis for expecting cloudiness to increase with inversion strength is well established: stronger inversions limit entrainment of warm and dry free-tropospheric air, keeping the BL shallower, cooler, moister, and hence cloudier. That CMIP3 GCMs do not respond this way (as evidenced by Fig. 5a and the fact that GCM cloud decreases as global warming increases EIS) suggests that subtropical low-cloud feedback may be systematically overpredicted in CMIP3 models.
Figure 9 shows dMLMcloud/dEIS computed from current-climate seasonal and geographic variability plotted against dMLMcloud/dEIS computed from differences between current- and future-climate mean values. In most models, projecting future cloud change from current-climate variability results in overpredicted cloud change. This suggests that other processes (e.g., direct CO2 effects, subsidence, or advective changes) have cloud-change damping impacts that do not project onto EIS change. Understanding these processes is important future work.
Correlations with other forcings are weaker and likely related to interconnection between forcings (as discussed in section 4). Divergence in particular does not seem to be a good predictor of intermodel differences in cloudiness or cloud change. This is in contrast to the predictions of CB09 and is probably due to the fact that 850-hPa-mean divergence changes tend to be weak in the CMIP3 archive as discussed above.
Changes in important MLM variables are presented in Table 6. Cloud base and BL depth (cloud top) changes are modest and vary in sign between regions. Excluding Australia, their combined effect is to increase cloud depth and hence liquid water path. This is consistent with the fact that we decreases robustly in all regions (probably due to EIS increase), moistening and cooling the BL. Even in the Australian region where cloud thickness decreases on average, mean LWP increases slightly due to stronger adiabatic liquid water increase with height (as noted by Somerville and Remer 1984).
6. Stratocumulus depletion mechanisms
Further insight into why the MLM predicts an increase in stratocumulus incidence under climate change is achieved by examining the conditions under which Sc is diagnosed as absent and how these conditions are affected by climate change. Figure 10a documents how frequently each no-Sc condition is triggered in the current climate. The main source of reduced cloudiness in the current climate is decoupling due to BL depth >2 km. Adding decoupling due to large-scale upward motion (D < 0) increases the frequency of decoupling by an additional ~20%. Shallow, clear-sky BLs driven by shear occur less than 10% of the time in all regions, with nonconvergent solutions occurring less than 5% of the time. Spread between model runs is largest for frequency of deep decoupled BLs, suggesting that this is the main source of intermodel and interregion spread.
Global warming–induced changes in Sc-free frequency (Fig. 10b) are more complicated. Changes in overly deep BLs still tend to be most important, with increased Sc incidence due primarily to BL depths exceeding 2 km less frequently. However, other mechanisms play an important role in some regions. Nonconverged cases are again negligible except in Canary Islands and Australian regions, where they cancel increases in shear-driven BLs and cases of upward motion. The frequency of BL depth >2 km decreases robustly in all regions, consistent with reduced entrainment as EIS rises. Changes in the frequency of clear, shear-driven BLs are regionally dependent, but most models show responses of the same sign within a given region. This suggests that the MLM may be capturing real regional differences, though the reason for this response is still unclear.
When large-scale conditions from CMIP3 models are used to drive a MLM, low cloud changes (calculated as the frequency of cloudy equilibrium solutions) act as a negative feedback on global warming. The primary mechanism for cloud change is a robust increase of ~20% in estimated inversion strength (EIS) in the driving models. The physical explanation is that, as inversion strength increases, entrainment of warm dry free-tropospheric air decreases, which causes liquid water path (LWP) to rise and, because the BL does not grow as deep, decreases the frequency of decoupled conditions. In contrast to previous studies (e.g., CB09), weakening divergence does not play a leading role (perhaps because mean divergence changes at 850 hPa are fairly weak in the GCMs studied).
CMIP3 GCMs generally predict positive cloud feedback—the opposite of that predicted by the MLM. The reason for this discrepancy is that GCM cloud fraction is not sensitive enough to EIS—both for determining climate change response and for predicting current-climate spatial and temporal variations. Since EIS is observed to have a strong positive correlation with present-day low cloudiness and low cloudiness has a large effect on surface temperature, our work suggests that global warming may be weaker than expected from the CMIP3 models. Boundary layer depth rise and apparent increase in BL stratification suggest that GCM behavior may result from the use of boundary layer schemes that are driven from the surface layer and, hence, do not depend properly on inversion properties. Because more models with explicit treatment of cloud-top cooling (e.g., Bretherton and Park 2009) are expected in the next CMIP intercomparison, we are excited to explore this with new data.
That these GCMs are not responding properly to thermodynamic forcing changes implies that CMIP3 problems arise from discrepancies in the physics parameterizations rather than dynamics [supporting the findings of Bony and Dufresne (2005) and Medeiros et al. (2008)]. Thus, one might expect a smaller spread in climate change predictions for our simulations, which are roughly analogous to runs from the various GCMs with their physics parameterizations replaced by the MLM. That intermodel spread in cloud fraction is not reduced in our simulations implies that improving physics parameterizations is a necessary but not sufficient condition for decreasing intermodel low-cloud spread. Spread in EIS in CMIP3 models is significant but does not influence the GCM Sc-region cloudiness changes because the included models are not sensitive enough to EIS changes. On this point, it should be noted that the UKMO HadGEM1 model, which Clement et al. (2009) found to best simulate current-climate relationships between cloud and inversion strength (among other things), is not included here because it lacks the needed MLM forcings. Interestingly, despite responding better to current-climate EIS changes, this model predicts positive low cloud feedback.
This paper points to the need for work on several fronts. Clearly, GCM representation of subtropical low clouds must be improved if the magnitude of future warming is to be accurately determined. Physical parameterizations in particular must better handle sensitivity to forcings. Determination of large-scale conditions must also be improved if model predictions are to converge. There is a good chance that next-generation GCMs are already better, so performing similar analysis on new data is a critical first step. Another important task is to understand why EIS increases in these models. If solid reasons for expecting EIS increase are uncovered, we will be able to say with a degree of certainty that subtropical low-cloud feedback will be negative.
We would like to acknowledge the modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the World Climate Research Program’s Working Group on Coupled Modeling for their roles in making available the CMIP3 multimodel dataset. Support of the CMIP3 dataset is provided by the U.S. Department of Energy (DOE) Office of Science. This work was performed under the auspices of DOE by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. All authors are supported by the DOE Office of Science’s Regional and Global Climate Modeling Program under the project “Identifying Robust Cloud Feedbacks in Observations and Models.” This study benefited from comments by Peter Blossey and one anonymous reviewer.
MLM Forcing Details
Many details must be considered in converting the limited daily CMIP3 data into mixed-layer model forcings. These are described below. Where possible, approximations are tested against exact solutions using California-region ERA-Interim data between 1990 and 2000. Correlations are calculated by considering each cell in the region and each time as independent.
Because only monthly averaged SST is available in the CMIP3 archive, we use the method of Taylor et al. (2000) to derive piecewise-linear daily SST, which reproduces the given monthly averages. Another missing variable is divergence, which we calculate as the mean over the layer between the surface and 850 mb. Divergence at each level is calculated from winds following appendix C of Holton (1992),
with eastward wind u, northward wind υ, and latitude φ and longitude λ (in radians). This calculation is inexact because the available output is coarser than used internally by the model. We test the induced error by comparing D output directly from ERA-Interim against that calculated from (A1). The two fields are correlated at r2 = 0.86 and their means differ by less than 5%. Variance, however, is reduced by 40%, probably owing to smoothing from using centered differencing for discretization.
Because the MLM uses a geometric height coordinate but CMIP data is only available on pressure levels, conversion between height and pressure levels is needed. We do this by calculating the 850-mb geopotential height Z850 and specifying vertical gradients in geometric height units. Z850 is computed using
where Rd is the gas constant for dry air (taken to be 287 J kg−1 K−1), g = 9.8 is the gravitational constant, and Tυ = T(1 + 0.61qυ) is the virtual temperature for temperature T and water vapor mixing ratio qυ. Applying our approximation to ERA-Interim data (for which Z850 is an output) yields negligible mean, variance, and correlation errors (<1%).
Computing free-tropospheric boundary conditions is another challenge. While it is easy to extract quantities at archived pressure levels, if the reference level is too low, it may be in the BL some of the time and, if it is too high, it may not be representative of air that would actually be entrained. To avoid this we construct synthetic “free-tropospheric 850-mb” values. We start by diagnosing the lowest free-tropospheric layer as the highest layer in the region where RH decreases rapidly with height, which we define as the highest layer with second derivative satisfying
where A and B are the mean and standard deviation of ∂2RH/∂p2 over the region below 600 hPa and α is a tunable parameter describing how distinct the base of the free troposphere must be from the background. This equation specifies how much stronger than background the curvature of the RH profile must be to still be in the inversion. Allowing the base of the free troposphere to move up is important for well-resolved profiles where several levels may have similar slope but is generally not active for coarse profiles. Since (A3) adds a thresholding effect, small changes in α can have large effects on individual profiles. In almost all cases, however, modest changes in α have no effect on diagnosed the free-troposphere base height. We take α = 0.2 because it gives reasonable results at both fine- and coarse-resolution for the first 90 days of 1990 using ERA-Interim data at 30°N, 130.5°W. Our method is illustrated in Fig. A1, and predicted free-troposphere bottoms are included as a white line in Fig. 2. We search for free-tropospheric layers below 600 mb. In some cases no inversion can be diagnosed. We consider such points to be D < 0 cases. Since lack of a BL is generally indicative of convective conditions, this is a reasonable approach. Other methods were explored and gave similar (but slightly worse) results.
If the diagnosed layer is higher than 850 hPa, we extrapolate a new 850-hPa value using least squares linear regression on the six lowest free-tropospheric layers. Values at 850 hPa are used directly if this level is not in the BL. Free-tropospheric slopes for the MLM are taken from the previously described least squares regression if 850 mb is in the BL; otherwise, slopes come from least squares regression on the six levels starting at 850 mb. The regression slope is then converted to height coordinates using the hydrostatic relation.
Horizontal advective tendencies in total water mixing ratio qt, T, and inversion height zi must also be derived. Because BL properties are closely tied to SST, we set ∇T = ∇SST and ∇qt = 0.8∇qs(SST, psurf), where qs is the saturation mixing ratio and Psurf is the surface pressure. The factor 0.8 accounts for the fact that oceanic near-surface RH is typically 80%. Use of daily averaged data to compute advection likely induces some error. Its magnitude cannot be ascertained because advection is not an ERA-Interim output. Computing ∇zi is more difficult. The lowest free-tropospheric layer computed above gives gradients in decent agreement with the ERA-Interim BL height output when all of its considerably larger number of vertical levels are used. When only CMIP3 levels are used, however, horizontal gradients end up being zero most of the time with occasional random jumps between levels. We also explored using the simple model of CB09 to make “first guess” snapshots of BL inversion height, but this also performed poorly. In the end, we concluded that information sufficient to compute BL gradients is not present in the CMIP3 archive, so we set BL advection equal to a constant value 0.49 mm s−1. This value comes from three months of satellite data for the southeast Pacific analyzed by Wood and Bretherton (2004) but is very similar to ERA-Interim data, which has a mean value of 0.35 mm s−1 when averaged over 1990–2000 and across all our regions except Australia (which is omitted because its height advection causes BL deepening). Rerunning our simulations with this seasonally and regionally varying BL height advection is shown in appendix B to make very little difference. For this reason, we do not think using constant BL height advection is a significant source of error in our study.
A drawback to modeling studies is that it is often unclear whether conclusions reflect underlying physical truth or arise from imperfections in the parameterizations. Testing model sensitivity to parameterization details is a good way to increase confidence in model results. To this end we have performed additional sets of simulations, summarized in Fig. B1. This figure shows the time-, space-, and model ensemble-averaged relationship between current-climate cloud fraction and EIS for a variety of sensitivity studies.
The impact of using seasonal and regionally varying BL height advection values from ERA-Interim is seen to be extremely small (red line). If boundary layer incursions above 850 hPa were not corrected (brown line), cloud fraction would increase slightly, but the functional relationship between EIS and cloud fraction would remain unchanged. Changing the limiting value for radiative cooling also has little impact on our simulations (Fig. 7), suggesting that our results are insensitive to our choice of radiative scheme.
Varying the depth beyond which decoupling occurs has a larger effect. Using a cutoff of 1 km [as suggested from the observational study of Wood and Bretherton(2004)] results in almost total depletion of cloud because our BL is usually deeper than 1 km. At a cutoff of 1.5 km, little cloudiness is found for EIS < 5 K, but the default slope is regained for stronger inversions. Increasing the cutoff to 2 (the default) or 2.5 km translates cloud fraction upward but maintains approximately the same slope. This suggests that our method’s ability to replicate observed sensitivity to EIS is independent of our choice of cutoff wherever matching the observed slope does not require negative cloud fraction.
Also included is a simulation with drizzle turned off. Because drizzle increases as clouds thicken and reduces BL growth through latent heating of the cloud layer, neglecting drizzle results in deeper BLs, which are more likely to be decoupled in our model. Although not included in our model, drizzle can also stratify the BL by heating the cloud layer and cooling the subcloud layer through evaporation (Caldwell et al. 2005). Including a more sophisticated decoupling criterion and allowing subcloud evaporation of drizzle may alter the impact of drizzle on our simulations in a way that we are not testing. Exploring this is important future work. Cloud response to EIS increase is damped in the nondrizzling run, which changes the functional relationship between cloud and EIS. To ensure our results are not dependent on our drizzle formulation, we repeated our climate change experiments with drizzle turned off (Fig. B2). Cloud responses with and without drizzle are very well correlated (r2 = 0.8), implying that the factors determining sensitivity are similar for both model formulations. As in the base case, the nonprecipitating case predicts generally more low cloud in the future. A difference between the cases is that intermodel spread is substantially smaller for the nondrizzly case (indicated by smaller spread in the y direction than in the x direction). This is unsurprising since drizzle was shown above to amplify cloud fraction response to EIS perturbations.
560 hPa was chosen to avoid ISCCP cloud height misplacement under strong inversions. Using 680 hPa gave similar results.
Daily anomalies from the climatological seasonal cycle show similar correlations.