A new coupled cloud physics–radiation parameterization of the bulk optical properties of ice clouds is presented. The parameterization is consistent with assumptions in the cloud physics scheme regarding particle size distributions (PSDs) and mass–dimensional relationships. The parameterization is based on a weighted ice crystal habit mixture model, and its bulk optical properties are parameterized as simple functions of wavelength and ice water content (IWC). This approach directly couples IWC to the bulk optical properties, negating the need for diagnosed variables, such as the ice crystal effective dimension. The parameterization is implemented into the Met Office Unified Model Global Atmosphere 5.0 (GA5) configuration. The GA5 configuration is used to simulate the annual 20-yr shortwave (SW) and longwave (LW) fluxes at the top of the atmosphere (TOA), as well as the temperature structure of the atmosphere, under various microphysical assumptions. The coupled parameterization is directly compared against the current operational radiation parameterization, while maintaining the same cloud physics assumptions. In this experiment, the impacts of the two parameterizations on the SW and LW radiative effects at TOA are also investigated and compared against observations. The 20-yr simulations are compared against the latest observations of the atmospheric temperature and radiative fluxes at TOA. The comparisons demonstrate that the choice of PSD and the assumed ice crystal shape distribution are as important as each other. Moreover, the consistent radiation parameterization removes a long-standing tropical troposphere cold temperature bias but slightly warms the southern midlatitudes by about 0.5 K.
The Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC; Stocker et al. 2013) concluded that the radiative coupling between clouds and the atmosphere remains one of the largest sources of uncertainty in understanding intermodel differences in predicting the equilibrium state of Earth’s climate. Stocker et al. (2013) concluded that most of the intermodel differences were due to low clouds. It is well known, however, that general circulation models (GCM) generally underpredict ice mass, compared to observations, by several or more factors (Waliser et al. 2009; Wu et al. 2009; Delanoë et al. 2011; Field et al. 2014). This general underprediction of high-troposphere ice mass by climate models creates problems when assessing the importance of different high cloud types to the climate sensitivity of GCMs. This is because the role of high clouds in determining the shortwave (SW) and longwave (LW) radiative exchanges may well be underestimated by those models (Baran 2012). Moreover, the parameterization of high clouds can affect the amount of low clouds predicted by the model through the vertical profile of radiative heating in the model, as shown by McFarquhar et al. (2003). To quantify the role of the ice cloud in the radiative coupling between the atmosphere and cloud, it is of primary importance to construct accurate parameterizations of its bulk optical properties. Unfortunately, this is currently far from being achieved, as different ice crystal models and parameterizations can lead to significantly different GCM-simulated tropical LW and SW flux differences at the top of the atmosphere (TOA) of between about 10 and 30 W m−2 and −10 and −35 W m−2, respectively (Baran 2012). The study by Fu (2007) showed that changing the aspect ratio of hexagonal ice columns from 1.0 to 0.1 leads to a shortwave cloud radiative effect (CRE) difference of about −40 W m−2 at TOA. The Fu (2007) result was obtained by assuming a completely overcast sky and an optical depth of 4.0. In general, a recent study by Williams and Webb (2009) showed that the divergence between different climate models in predicting the CRE of cirrus can be ±20 W m−2 and between −50 and 15 W m−2, in the tropics and extratropics, respectively.
However, the differences highlighted above could be due to different cloud physics parameterizations, such as particle size distributions (PSDs) and/or ice optics, being assumed in different climate models. The study by Edwards et al. (2007, hereafter E07) showed that it is possible to obtain insignificant differences in the simulated SW and LW fluxes at TOA between seemingly different ice crystal optical parameterizations if the PSD and shape of ice crystals are known and the same assumptions are made between the ice crystal size and environmental temperature. One further advance on this work would be to assume the same mass–dimensional relationship and PSD in the cloud physics and radiation schemes of the climate model. A second advance on the previous work would be to insist that the ice crystal model assumed in the radiation scheme must follow the currently observed area and mass–dimensional relationships that have been derived from in situ measurements (Mitchell et al. 2008; Baran 2009, 2012; Cotton et al. 2013).
The studies of (Baran and Labonnote 2007; Mitchell et al. 2008; Baran et al. 2009, 2011a,b, 2014a) have shown if the idealized ice crystal model follows observed area and mass–dimensional relationships, then the need for assumptions about relationships between an effective dimension of the ice crystal population (the effective dimension is defined below) and environmental temperature can be negated when predicting the radiative properties of ice cloud. The approach of the previous authors means that microphysical assumptions in both the cloud physics and radiation schemes in the climate model are consistent, whereas optical property parameterizations that rely on assumed relationships between an effective dimension and environmental temperature tend to be physically inconsistent between the cloud physics and radiation schemes (Baran 2012). This physical inconsistency arises because different assumptions are made between the cloud physics and radiation schemes with regard to shape of the PSD and the area and mass–dimensional relationships. In general, current parameterizations of ice crystal optical properties tend to be functions of ice crystal effective dimension, ice water content (IWC), and/or environmental temperature.
The ice crystal effective dimension De is defined as the ratio between the IWC and the volume extinction coefficient, where the volume extinction coefficient is the sum of the scattering and absorption coefficients (Francis et al. 1994, 1995). However, there are multiple definitions of De that can be employed, as shown by McFarquhar and Heymsfield (1998), and these definitions can vary in their values by much more than 20%. For some definitions of De, it is itself an optical property that can be used to describe the extinction of light in a cloud. This is because it is inversely related to the mass extinction coefficient (Mitchell et al. 1996; McFarquhar and Heymsfield 1998; Wyser and Yang 1998; Kokhanovsky 2004; Mitchell 2002). This inverse relationship is only valid for solar wavelengths, because at those wavelengths, where the wavelength of the incident solar radiation is much smaller than the largest dimension of the ice crystal, the extinction coefficient is twice the ice crystal orientation-averaged projected area (van de Hulst 1957). Therefore, given De and IWC, the radiative properties of cirrus can be uniquely determined for solar wavelengths under the assumption of a homogeneous distribution of IWC within a single layer.
Most of the current parameterizations of cirrus bulk optical properties used in climate models depend on empirically derived deterministic relationships between De and environmental temperature and/or IWC (Fu et al. 1999; Kristjánsson et al. 2000; E07; Bozzo et al. 2008; Hong et al. 2009; Gu et al. 2011; Yi et al. 2013). Generally, there is no evidence for exact deterministic relationships between De and environmental variables and/or IWC. In the studies reported by Lynch et al. (2002) there appear to be very weak correlations between De and environmental variables, such as in-cloud temperature. A further limitation of some of the earlier parameterizations cited above and the lack of correlation between De and in-cloud temperature and/or IWC, is that the historical in situ–measured PSDs on which they depend are now known to be biased toward populations of small ice crystals. This is due to the shattering of ice crystals on the housing and/or at the inlets of the microphysical probes that have been used to measure the size spectra (Field et al. 2003, 2006; Korolev et al. 2011, 2013). By mitigating the effects of shattering, more recent studies by Mitchell et al. (2011b) have demonstrated stronger correlations between De, in-cloud temperature, and IWC.
How the issue of shattering affects previous cirrus radiative parameterizations depends on the definition of De used in those parameterizations. For instance, it can be seen in Korolev et al. (2013) that shattering on microphysical probes affects, for particles of size less than about 500 μm and in the temperature range from 0° to −35°C, both the extinction and IWC in an approximately similar way (this is because extinction and IWC may have a similar dependence on the diameter of the ice crystals). The effect of shattering on the extinction and IWC is to increase them by, on average, 20% and 30%, respectively (Korolev et al. 2013). Therefore, the effect of shattering on their ratio will be small but may not completely cancel the effect on estimates of De. However, if De is defined as the ratio of the third moment (i.e., volume) to the second moment (i.e., area) of the PSD, then shattering will bias the calculated value of De to smaller values. The latter definition of De is the one that is assumed in the current Met Office operational Global Atmosphere model, and the PSDs used in the current Met Office parameterization were applied without removing any shattered ice crystal artifacts. In the study by Mitchell et al. (2010), it is shown that shattering could cause the in situ–estimated De to be underpredicted by several factors, and in that paper for ice crystal size less than 240 μm, the mass of ice tends to that of a sphere. The impact of shattering on IWC and extinction is still uncertain, as it does depend on the extent of the PSD and on what size of ice crystal contributes most to the IWC. Clearly, further research in this area is required to determine the small ice mode (Korolev et al. 2013). A possible consequence of this underprediction of De is to bias the averaged mass extinction of clouds to higher values, which means that the reflected SW radiation at TOA is biased to larger values (since mass extinction is inversely related to De). Moreover, in the SW, the underprediction of IWC in most GCMs leads to darker clouds, while the bias toward small De leads to brighter clouds. Hence, the two errors effectively cancel, which may lead to about the right emergent SW fluxes being predicted in the GCM. This cancellation of SW error excludes the possibility of exposing systematic climate model errors by comparing model SW flux predictions directly against SW flux measurements. However, there are large uncertainties in a global climate model, so this SW effect may not be apparent. It is more difficult to state the impact of the parameterization error on the outgoing LW flux, as that fundamentally depends on where the model locates the cloud in the atmosphere. For instance, if the model predicts the cirrus to be too optically thin and not high enough in the atmosphere, then this would increase the outgoing LW flux, and the overestimate of mass extinction may not necessarily be sufficient to cancel the error. To assess the performance of models in the SW and LW, it may prove more useful to move away from flux comparisons to high-resolution radiance comparisons, as suggested by Goody et al. (1998). However, to test model physical consistency across the electromagnetic spectrum, SW and LW high-resolution radiance measurements need to be combined (Baran and Francis 2004).
Yet a further difficulty with previous parameterizations is that, generally, the assumed PSDs used in the cloud physics and radiation schemes are not consistent. This means that De-based radiative parameterizations result in different De value being predicted in the radiation scheme, while in the cloud physics scheme, for the same mass of ice, a different De could be predicted. If so, this is physically inconsistent and, as such, should be removed from climate models. Note that the inconsistent types of parameterizations discussed above are prevalent in climate models that are used in the IPCC series of reports.
For De to be generally useful, it must be applicable across the spectrum. That is, the radiative properties of ice cloud must be uniquely determined from De and IWC alone. Unfortunately, this is not generally found to be the case. At infrared and radar wavelengths, the concept of De is no longer valid, as demonstrated and discussed by several authors (Mitchell 2002; Baran 2005, 2007; De Leon and Haigh 2007; Mitchell et al. 2011a; Baran et al. 2011b). The reason for the breakdown in the validity of the De concept at infrared wavelengths is simply because, at those wavelengths, the ice crystal size becomes small relative to the incident wavelength, which means that the mass extinction coefficient is no longer inversely related to De (Baran 2005). More fundamentally, as the ice crystals become small relative to the incident wavelength at infrared wavelengths, the physical process of absorption through electromagnetic wave resonance effects becomes more important, so absorption by ice is no longer uniquely determined by the ice crystal volume or projected area (Mitchell 2002; Mitchell et al. 2011a). At radar wavelengths, the radar reflectivity becomes more sensitive to aggregated ice crystals. As new ice crystals grow in size by vapor diffusion and also become aggregated, with the larger ice crystals more likely to aggregate, the ice crystal mass–dimensional relationship asymptotes to D2, where D is the ice crystal maximum dimension (Westbrook et al. 2004; Heymsfield et al. 2010; Cotton et al. 2013). As a result of this asymptotic behavior in ice mass, the dependence of the numerator in the De formula tends to D2, while the denominator is also proportional to D2; hence De in regions of ice crystal aggregation, tends to constant values. This behavior in De can be seen in some of the profiled aircraft observations of Francis (1995), which show little vertical variation in effective radius (i.e., effective radius = De/2), and it is noted by Francis (1995) that the lack of variation in effective radius was due to the similar dependence that mass and area have on ice crystal size. This same behavior was also noted by McFarquhar and Heymsfield (1998). Clearly, in this situation, De tends to a constant value with respect to the ratio of IWC to extinction. Moreover, because of the nonuniqueness of cirrus PSDs, there is the possibility that there is no unique radiative solution for a particular value of De (Mitchell 2002).
In this paper, an alternative to the De approach is presented for the parameterization of the bulk optical properties of ice cloud. That is, the bulk optical properties are directly linked to the climate model prediction of IWC. This alternative approach enables the radiative properties of ice clouds to be directly related to a GCM prognostic variable rather than some generally diagnosed quantity, such as De. A similar approach was used by McFarquhar et al. (2003), where the bulk optical properties were equivalently parameterized as a function of IWC.
The climate model used as a basis for the simulations in this paper is the Met Office Unified Model Global Atmosphere 5.0 (GA5) configuration. The new parameterization is tested by comparing the results of the climate model simulations against the Clouds and the Earth’s Radiant Energy System (CERES) reanalyzed product of the annual 20-yr globally averaged reflected and outgoing longwave (OLR) fluxes at TOA. The CERES product used here is described in Loeb et al. (2009). The parameterization proposed by E07, a De-based scheme, is also compared against the same CERES products under the same microphysical assumptions as used in the experiments using the new parameterization. The impacts of all radiation parameterizations on the temperature structure of the troposphere are also discussed and are compared against the Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim) of atmospheric temperature observations, covering the period 1989–2008 (Dee et al. 2011).
The paper is split into the following sections. Section 2 briefly describes the GA5 configuration, and the atmospheric parameterizations within it that are most pertinent to this paper. Section 3 briefly describes the self-consistent scattering model for cirrus and describes the bulk optical properties that the model predicts. Section 4 presents the parameterization along with an error analysis of the parameterization, and section 5 describes the results of the various experiments. Section 6 presents the conclusions and discusses the most important points of this paper.
2. The GA5 configuration and atmospheric parameterizations
The approach of the Met Office in simulating the weather and climate is to keep the atmospheric physics parameterizations used for both the same, as far as possible, for use on all time scales (Brown et al. 2012). This has led to an annual definition of a Global Atmosphere (GA) configuration of the Met Office Unified Model, and this definition of a single GA model has been described by Walters et al. (2011). In this paper, use is made of the GA5 configuration, which is based on the GA4 configuration (Walters et al. 2013) but includes a new dynamical core. The new dynamical core is fully described by Wood et al. (2014). The details of the atmospheric parameterizations used in the GA5 configuration can be generally found in Walters et al. (2013). The GA5 configuration used in this study consists of 85 levels, the highest level being at an altitude of 85 km, and all the simulations used in this study are run at a horizontal resolution of about 135 km in the midlatitudes.
However, the GA parameterizations that are particularly pertinent to this paper are briefly discussed here. The radiative transfer model used in the GA series of models is the Edwards and Slingo (1996) scheme. Gaseous absorption is treated within the radiative transfer model by using the correlated-k method with six bands in the shortwave and nine bands in the longwave. In the longwave and shortwave, multiple scattering is fully included in the Edwards and Slingo (1996) scheme. A further difference between GA4 and GA5 is that, in GA4, a full radiative transfer calculation is done every three hours, and a fast calculation is done every hour using a few k terms to account for cloud changes, as described in Manners et al. (2009), whereas, in the GA5 configuration, a full radiative transfer calculation is done every hour. For a more detailed description of the radiation scheme, see section 2.2 of Walters et al. (2013), and the cloud scheme is fully described in sections 2.3–2.4 and in 3.3–3.4 of that paper. The GA convection scheme is described in section 3.6 of Walters et al. (2013).
The current GA operational ice optics parameterization is described by E07. The bulk ice optical parameterization developed in that paper is based on the 8-monomer hexagonal ice crystal model developed by Yang and Liou (1998). The hexagonal ice aggregate model is invariant with respect to size, so the asymmetry parameter will also remain invariant with respect to size. Observations generally show that aggregates of ice crystals grow longitudinally rather than radially, so their area ratios (i.e., the ratio of the projected area of a nonspherical particle to the area of the circumscribing circle of the same maximum dimension as the nonspherical particle) decrease as a function of increasing ice crystal size (Heymsfield and Miloshevich 2003; Field et al. 2008; McFarquhar et al. 2013). The asymmetry parameter is naturally related to the ice crystal size; so it is also related to their area ratio or aspect ratio (Fu 2007); therefore, this parameter should change with ice crystal size. Furthermore, observations used by Baum et al. (2005) and Baran et al. (2011b) show that the hexagonal ice aggregate model overpredicts in situ–derived IWC by several factors. The hexagonal ice aggregate shortwave and longwave bulk optical properties used in the E07 parameterization were compiled by Baran and Francis (2004). To address the physical inconsistencies inherent in the hexagonal ice aggregate model discussed above, an ensemble model of ice crystals developed by Baran and Labonnote (2007) has been proposed, and this model is briefly described in section 3 of this paper.
The assumed ice PSD that is currently used in GA5 is a parameterization developed by Houze et al. (1979); it is based on 37 in situ PSDs, and 90% of these were measured at temperatures warmer than −30°C. Currently, the Houze et al. (1979) estimated PSD is kept constant at temperatures colder than −35°C. This assumption means that at much colder temperatures, there will be orders of magnitude more frequently occurring large ice crystals than there should be. This has clear implications for the assumed fall speed, which must be artificially increased to accommodate space-based radiometric SW and OLR measurements.
A more recent parameterization of the ice PSD has been proposed by Field et al. (2007, hereafter F07). The parameterization of F07 is a moment estimation scheme, and it is based on 10 000 in situ PSDs that were measured between the temperatures of about 0° and −60°C in the midlatitudes and tropics. Figure 1 shows the measurement space in IWC and cloud temperature over which the F07 and Houze et al. (1979) parameterizations are based. The figure shows that the F07 PSD parameterization is based on measurements over a far greater range of IWC and cloud temperature than the Houze et al. (1979) parameterization. As previously stated, the historical PSDs have been affected by ice crystal shattering. The technique of filtering employed by F07 to reduce the effect of shattering is described in Field et al. (2006) and so will not be repeated here. The F07 parameterization also ignores all in situ measurements of ice crystals with maximum dimensions less than 100 μm. For ice crystals less than 100 μm in size, the F07 parameterization assumes an exponential distribution (Field and Heymsfield 2003), which is added to a modified gamma distribution at an ice crystal size of about 100 μm. As previously discussed, shattered ice crystal artifacts may still be present for ice crystals less than 500 μm in size (Korolev et al. 2011, 2013). However, the extent to which the F07 parameterization has been affected by shattering is unclear. In Korolev et al. (2013) it is shown that if filtering alone is applied to remove shattered artifacts (Field et al. 2003, 2006), then significant divergence from the best estimate of the PSD does not occur until ice crystal size is less than about 200–175 μm. It cannot be categorically stated that the F07 parameterization has not been affected by shattering, and if so, by how much, since the parameterization is based on in situ data that covers a far greater temperature and altitude range than was considered in Korolev et al. (2013). However, the F07 parameterization is still a more representative parameterization of the PSD than Houze et al. (1979), and the latter parameterization did not employ any techniques to try to remove shattered ice crystal artifacts from their parameterization. Moreover, Furtado et al. (2014) show that the F07 parameterization is a better fit to in situ–derived higher moments of the PSD than the Houze et al. (1979) parameterization. The in situ data used by Furtado et al. (2014) is based on microphysical probes, which have the necessary modifications attached to them to reduce shattering, and filtering has been applied to further reduce the effects of shattering on the moments of the PSD (Field et al. 2003, 2006; Cotton et al. 2013; Korolev et al. 2011, 2013).
The moments in the F07 parameterization are predicted as a function of IWC and cloud temperature, where the moment carrying ice mass is generally related to the IWC, and all other moments are generated from the second moment and cloud temperature. To generate the PSD from the F07 moment estimation parameterization, the second moment is related to an assumed mass–dimensional relationship.
In GA5, the Brown and Francis (1995) derived mass–dimensional relationship is assumed, where the mass of each ice crystal is given by 0.0185D1.9, where D is the ice crystal maximum dimension (both mass and D are in SI units, which are used throughout this paper). However, a recent midlatitude cirrus campaign conducted by Cotton et al. 2013 found that the Brown and Francis (1995) relationship overestimated in situ bulk measurements of IWC by, at most, a factor of 2. In the same paper, the authors found that the mass–dimensional relationship that best fitted the measurement of bulk IWC was found to be 0.0257D2. Moreover, the mass–dimensional relationship derived by Cotton et al. 2013 is also consistent with the mass–dimensional relationship found independently by Heymsfield et al. (2010) using bulk measurements of IWC obtained in midlatitude and tropical cirrus. The ice microphysics assumed in the model experiments used in this paper are based on the F07 PSD parameterization, the mass–dimensional relationship derived by Cotton et al. 2013, and a fall speed parameterization that has been developed by Furtado et al. (2014). In the next section, the ensemble model of Baran and Labonnote (2007) is described, along with its bulk optical properties.
3. The ensemble model and its single-scattering properties
It has recently been argued by Baran (2009, 2012) that ice optical properties should be directly linked to global model prognostic variables, such as IWC. It was also argued by the same author that the same PSDs should be assumed in both the radiation and cloud physics schemes within climate models. To this end, the ensemble model of cirrus ice crystals was developed to predict the orientation-averaged cross sections of naturally occurring ice crystals and conserve ice crystal mass (Baran and Labonnote 2007). The ensemble model of cirrus ice crystals is shown in Fig. 2, and the figure shows that the ensemble model consists of six elements, the first of which is the hexagonal ice column of aspect ratio unity; the second element is the six-branched bullet rosette. Thereafter; the ensemble model consists of hexagonal ice aggregates forming 3-, 5-, 8-, and 10-monomer ice aggregates. Each of the ice aggregates is constructed by arbitrarily attaching hexagonal monomers to each other until an ice aggregate is constructed. The individual monomers are sufficiently spaced from each other to minimize multiple reflections between monomers. The first element represents the smaller sizes of ice crystals in the PSD, while the hexagonal ice aggregates represent the process of ice crystal aggregation and, thus, represent the larger sizes of ice crystals in the PSD. The members of the ensemble are distributed into six equal size intervals of the F07 PSD. This distribution means that, for radiative transfer in the solar and thermal regions of the spectrum, the cloud extinction is chiefly determined by the first element, with no contributions from the fourth member onward; however, for microwave and radar calculations, the cloud extinction would be weighted toward the more aggregated ensemble members (i.e., fourth member onward). For the same ice crystal maximum dimension, the more aggregated ensemble members would have much lower values of orientation-averaged cross sections than the first member of the ensemble model. The weightings critically depend on the shape of the PSD, and the weightings could change, given different shapes of the PSD. The impact of changing the weighting of the ensemble model on the reflected shortwave and OLR at TOA and on the temperature structure of the troposphere, is explored in one of the experiments contained in this study.
With the distribution of ensemble members described above, a series of papers by Baran et al. (2009, 2011a,b, 2014a) demonstrated that the model could predict, to within current experimental uncertainties, the solar volume extinction coefficient and IWC of midlatitude and tropical cirrus. More recently Baran et al. (2014b) demonstrated that the ensemble model prediction of the area ratio was within the range found by a number of independent authors, who used in situ estimations of the area ratio obtained in midlatitude, tropical, and Arctic cirrus (McFarquhar et al. 2013; Field et al. 2008; Heymsfield and Miloshevich 2003).
Other habit mixture models have also been proposed, such as those of Baum et al. (2005, 2011), that incorporate other observed habits not considered here, such as hexagonal ice plates, aggregates of plates, hollow columns, and hollow bullet rosettes.
Moreover, the ensemble model has been demonstrated to be physically consistent across the electromagnetic spectrum, from the UV to the radar frequency of 35 GHz, using consistent microphysics across the spectrum (Baran et al. 2011b, 2014a). Using measurements from the UV-to-radar frequencies, it was shown by Baran et al. (2014a) that for one case of optically thin midlatitude cirrus (i.e., solar optical depth ≪1), the ensemble model could predict aircraft-mounted backscatter lidar-derived UV volume extinction coefficients to generally within ±25% for all altitudes considered. Moreover, it was shown that aircraft-based high-resolution brightness temperature measurements at wavelengths between about 8.0 and 12.0 μm and between about 3.3 and 5.0 μm could be simulated to within ±1 and ±2 K of the measured brightness temperatures, respectively. Therefore, the optical properties used in the parameterization presented in this paper have been demonstrated to be physically consistent across those parts of the electromagnetic spectrum of relevance to climate models. Furthermore, a recent global study by Vidot et al. (2014, manuscript submitted to J. Geophys. Res.), demonstrated that the ensemble model ice optical properties produced a mean departure of only 0.43 K from space-based brightness temperature measurements. The measurements were located in the terrestrial window region at 8.0, 11.0, and 12.0 μm. This study was comprised of 26 791vertical profiles of cirrus, retrieved from radar and lidar profiles, with visible optical depths ranging from 0.03 to 4.0, and these profiles were distributed between the latitudes of about 70°N and 60°S, and at altitudes from about 440 to 50 hPa. Therefore, given the range of cirrus used in the study by Vidot et al. (2014, manuscript submitted to J. Geophys. Res.), the ice optical properties used here can be applied to a climate model.
The consistent microphysics used in Baran et al. (2014a) was derived from the ensemble model, and the derivation is described in Baran et al. (2011b). In that paper, using observational data of ice crystal mass obtained in deep tropical convection by Field et al. (2008), the ensemble model was shown to fit a mass–dimensional relationship of the form 0.04D2. This mass–dimensional relationship was shown by Cotton et al. (2013) to be within the experimental uncertainty of other independently derived relationships, which were based on in situ bulk measurements of IWC obtained in a variety of cirrus. In this study, the ensemble model mass–dimensional relationship is used to generate the F07 PSDs, and the impacts of these PSDs on the reflected shortwave and OLR at TOA, as well as on the temperature structure of the troposphere, are discussed in a later section of the paper.
To simulate the SW and LW fluxes in the climate model the total optical properties are required, which are the bulk extinction coefficient βext, bulk scattering coefficient βsca, single-scattering albedo ω0, and the asymmetry parameter g. These are wavelength dependent, but this dependency has been dropped here for reasons of clarity. The methods used to compute the total optical properties predicted by the ensemble model have already been described by Baran et al. (2014a) and will not be repeated here for reasons of brevity. The bulk optical property βext/sca is defined by the following equation:
where Cext/sca() are the orientation-averaged scattering and extinction cross sections of each of the elements of the ensemble model, the vector l contains each ensemble model element as a function of its size, and n() is given by the F07 moment estimation parameterization. Therefore, the bulk extinction and scattering coefficients can be expressed as functions of IWC and temperature, since these two parameters are used to generate the F07 PSDs. An alternative to this, if one wanted to relate the bulk optical properties to some size dependency would be to obtain from the F07 PSDs the ratio between the third and second moments. This would then express the bulk optical properties as functions of the mean mass-weighted size of the PSD. However, here the intention is to relate directly the climate model prognostic second moment [i.e., IWC, because m(D) = constant × D2] to the bulk optical properties without the need to use diagnosed variables. The methodology adopted here, as previously mentioned, is equivalent to that of McFarquhar et al. (2003); in that paper, the effective radius of the PSD was expressed explicitly as a function of IWC, and then the bulk optical properties were themselves expressed as a function of the effective radius. This is equivalent to expressing the bulk optical properties as a function of IWC, which is the intention here.
Given βext and βsca, then ω0 can be found from
and the bulk asymmetry parameter g is defined as follows:
where all the terms have been previously defined.
Equations (1)–(3) were used to compute the ensemble model bulk optical properties at 145 wavelengths between 0.2 and 120 μm. To generate the F07 PSDs, the ensemble model mass–dimensional relationship was assumed. The PSDs themselves were estimated from a database consisting of 20 662 estimates of IWC and in-cloud temperature measurements obtained from a number of cirrus field campaigns, which were located in different regions of the world. The construction of the 20 662 PSD database is fully described by Baran et al. (2014a). The latest estimates by Korolev et al. (2013) suggest that the impact of shattering on in situ derivations of IWC using historical PSDs could increase those estimates by up to about 30% at the heights and temperatures they considered. This does mean that some of the IWCs used here to generate the F07 PSDs could be similarly biased. However, the previously discussed remote sensing studies using the ensemble model optical properties at wavelengths relevant to a climate model indicate that the forward modeled measurements are generally within the uncertainties of those active and passive measurements. The wavelength resolution at which the ensemble model bulk optical properties were calculated is shown in Figs. 3a,b. The figures show the real and imaginary indices of solid ice as compiled by Warren and Brandt (2008), and the circles superimposed on the refractive indices are the wavelengths at which the bulk optical properties were calculated. It can be seen from Figs. 3a,b that the wavelength resolution used in this paper is sufficient to capture the rapid changes in the ice refractive index, especially between 1.5 and 4.0 μm, and at wavelengths in the terrestrial window and far-infrared regions.
In the Met Office Unified Model, clouds are represented by vertical profiles of their mixing ratios with respect to air. Therefore, the volume extinction and scattering coefficients become the mass extinction and mass scattering coefficients. The spectral variation of the mass scattering coefficient Ksca and the asymmetry parameter g are shown as a function of qi (i.e., the ice mass mixing ratio) in Fig. 4a and Fig. 4b, respectively, where, in Fig. 4a, Ksca is a product of ω0 and Kext. The figure shows the variation of Ksca and g over the range of IWC used in the parameterization; this range goes from about 0.004 kg kg−1 to about 9.0 × 10−9 kg kg−1. These seven values of IWC were chosen from the possible 20 662 values to qualitatively illustrate the dependence of Ksca and g on the shape of the PSD as a function of qi. Figure 4a shows the spectral variation of Ksca, and it can be seen from the figure that at the highest values of qi, Ksca is only weakly dependent on wavelength; this is because, at the highest values of qi, the PSDs are very broad. Conversely, at the lowest values of qi, the spectral variation of Ksca is significant, because the PSDs at those values are narrow. The figure also shows the wavelengths at which absorption is highest, such as between about 2.0 and 3.5 μm and at 10.0 μm and then again at around 40 μm. Note, interestingly, in the far-infrared, Ksca increases between 20 and 30 μm because of the imaginary index of solid ice decreasing between those wavelengths. Figure 4b shows the spectral variation of g as a function of wavelength and IWC. At the lowest wavelength of 0.20 μm, the asymmetry parameter has a value close to 0.70, and at the wavelength of 100 μm, g can decrease to values less than about 0.2. This behavior is expected at the longer wavelengths, as the ratio between ice crystal size and wavelength is small, and the IWCs are also very low, meaning very narrow PSDs. At around 3.0 μm, the asymmetry parameter can increase to near 1.0; this is due to the very high ice absorption that occurs at that wavelength. To be computationally efficient in a climate or numerical weather prediction model, the simplest parameterizations of Kext, ω0, and g are sought, but as can be seen from Fig. 4a and Fig. 4b, the variation of the bulk optical properties may be difficult to capture by simple parameterizations. The parameterization of the bulk optical properties is discussed and its accuracy tested in the section that follows.
4. The parameterization of the bulk optical properties
To first order, the mass extinction and mass scattering coefficients are the most important terms to be parameterized accurately to simulate the radiative transfer of solar and infrared radiation through cirrus. It can be seen from Fig. 4a that Ksca is only weakly dependent on wavelength, at the highest values of qi, and that it only becomes strongly dependent on wavelength when qi < 10−4 kg kg−1. This dependence of the first-order terms on wavelength suggests that a simple linear parameterization of these terms could suffice. As previously mentioned, the database consists of 20 662 values of qi, and for each of these values, there are 145 bulk optical properties. The total number of points in the database is, therefore, 2.996 × 106. From this database, a total number of 28 values of qi were selected to capture the monotonically decreasing values of qi in near-regular spacing. The values of the qi 28 points ranged from 0.004 to 10−9 kg kg−1. With monotonically decreasing values of qi, simple fits were found for Kext, ω0, and g as a function of wavelength and qi. The parameterized bulk optical properties are described by the following equations:
where the form of the above equations was estimated by trial and error; the prefactors shown in Eqs. (4), (5), and (7) were found by nonlinear least squares fitting; and λ is the wavelength. The values estimated for each of the prefactors and exponents are listed in Table 1 and Table 2, which were generated assuming the mass–dimensional relationship derived by Baran et al. (2011b) for the six shortwave and nine longwave bands, respectively. In Table 1 and Table 2, the wavelengths in each of the bands were not regularly spaced, as can be seen in Fig. 3a and Fig. 3b, so the parameterizations were implicitly weighted toward the part of each band that had the most wavelengths. The bulk optical properties described by Eqs. (4)–(7) were not band-averaged by the top-of-atmosphere solar and terrestrial irradiances. This is because band averaging should not be necessary, as the model is well resolved in spectral space, and errors in the parameterization and GCM will be more significant than differences in the bulk optical properties between their nonaveraged and band-averaged values. However, this latter statement was tested by averaging the bulk optical properties over the solar and terrestrial irradiances using the methodology of Lindner and Li (2000). It was found that the averaging had a small effect on the prefactors and exponents presented in Tables 1 and 2 (not shown here for reasons of brevity), so its omission will not change any of the conclusions reached in this paper. The band limits shown in the tables are currently assumed in the GA configurations. It should be noted here that the parameterizations shown by the above equations are not fixed to the bands listed in the tables because they are functions of wavelength and can be applied to any of the GA configurations or any other climate model.
The ensemble model mass–dimensional relationship has been used here because it is consistent with the lidar and radiometric observations of semitransparent cirrus presented in Baran et al. (2014a). However, it is yet to be confirmed whether the mass–dimensional relationship used in Baran et al. (2014a) is applicable over a greater range of optical depth. The new bulk optical property parameterization is now compared against the E07 parameterization for three values of the IWC or qi and two values of the environmental temperature. The three IWC and two temperature values assumed are 1.0 × 10−3, 1.0 × 10−5, and 1.0 × 10−7 kg kg−1, and −60° and −30°C, respectively. As a reminder, the E07 parameterization assumes a relationship between De and the environmental temperature, and De is itself inversely related to the mass extinction coefficient. Therefore, at the coldest temperature, the mass extinction coefficient will be about twice as great as the mass extinction coefficient at the warmest temperature. These two temperatures will provide a range of mass extinction coefficients over which the new parameterization can be directly compared against that of E07. The range chosen for the ice mass mixing ratio is typical of what might be found in a climate model. The two bulk optical property parameterizations are compared against each other between the wavelength intervals of 2.0 × 10−7 to 3.2 × 10−7 m and 8.33 × 10−6 to 1.25 × 10−5 m, which are shortwave band 1 and longwave band 5 in the GCM. These two bands are chosen because shortwave spectral band 1 is largely nonabsorbing and longwave spectral band 5 covers the terrestrial window region as shown by Tables 1 and 2. The results of the comparisons are shown in Figs. 5a–c, where Fig. 5a shows the results for Kext for shortwave spectral band 1, while Figs. 5b,c show the results for Kabs (which is simply Kext − Ksca) and the asymmetry parameter for longwave spectral band 5, respectively. In the figures, the results obtained at the cold and warm temperatures using the E07 parameterization are shown by the open and filled circles, respectively, and the new parameterization is shown by the diagonal cross symbol. Figure 5a shows that the E07 parameterization-predicted mass extinction coefficient will, at the coldest and warmest temperatures, be greater than Eq. (4) by factors of almost 2 and just over 1, respectively. At shortwave spectral band 1, the single-scattering albedo is very close to 1 for both parameterizations, and the g value calculated using the E07 parameterization does not vary and is fixed, as a function of IWC, at a value of about 0.74. However, Eq. (7) varies between 0.77 and 0.74, between the highest and lowest values of IWC, respectively. This is because, at the highest value of IWC, the PSD is broader and will contain a greater occurrence of larger ice crystals than the smallest assumed value of IWC. This change in the shape of the PSD will naturally lead to larger and smaller g values, respectively. Because of the E07 parameterization being based on an invariant ice crystal model with respect to ice crystal size, at nonabsorbing wavelengths, its predicted g values cannot change irrespective of assumed temperature and/or IWC values.
Figure 5b shows that for the spectral band in the terrestrial window region, the E07 parameterization is about 50% more absorbing than the new parameterization at the lowest temperature at all values of IWC considered. However, at the temperature of −30°C, the new parameterization is about 10% more absorbing than the E07 parameterization. The change in ω0 between the two temperatures using the E07 parameterization is largely invariant with IWC, changing from 0.50 to about 0.52, respectively, while Eq. (6) remains invariant with IWC and predicts a value of about 0.54 for all values of IWC considered. Figure 5c shows the expected invariance in the asymmetry parameter calculated using the E07 parameterization. However, the asymmetry parameter calculated using Eq. (7) decreases by about 6% between the highest and lowest values of IWC. This change in the asymmetry parameter reflects the change in the shape of the PSD, as previously discussed. The comparisons between the two parameterizations tend to show, at least at cirrus forming temperatures, that the E07 parameterization, for the same ice mass at the colder temperatures, should increase and decrease the outgoing SW and OLR relative to the new parameterization, respectively.
The accuracy of each of the parameterizations shown in Eqs. (4)–(7) is now examined using the full database of 2.995 × 106 values, minus the 28 values that have already been used to obtain the parameterizations. The error is calculated as the relative percent error, defined by the following equation:
where the parameters trueλ and estimateλ are the exact bulk optical properties obtained at each wavelength from the 20 662 PSD database, which consists of 2.995 × 106 values and their estimated values using Eqs. (4)–(7), respectively. The dependence of ε(λ) on qi has been dropped from Eq. (8) for reasons of clarity. The results found for ε(λ) are presented as a function of λ and qi. The distribution of ε(λ) in qi and λ space is shown in Figs. 6a–c, respectively. The figures reflect, to some extent, the change in the imaginary refractive index of ice. At wavelengths in the near-infrared and window regions, at around 2.0–3.0 μm, and at wavelengths greater than about 30 μm, the ε(λ) errors in Kext and ω0 are greatest, which is where solid ice is most absorbing. However, the errors are generally greatest for qi < 10−6 kg kg−1. In general, in situ measurements of IWC down to such values are not, at the moment, possible, and current global model predictions of qi are only occasionally less than 10−6 kg kg−1. Moreover, the first-order term is Kext, which is observationally not known to within a conservative ±50% (Baran et al. 2009).
The relative errors in the parameterization for the asymmetry parameter at the longer wavelengths in the far-infrared generally follow Figs. 6a and 6b. However, in the important shortwave region and terrestrial window regions, the errors are within a few percent for all qi values considered. However, for λ ≥ 40 μm, the errors in the g parameterization generally increase, with the largest errors occurring when qi < 10−6 kg kg−1. The relative errors found for g in the shortwave will not significantly bias the global model results. For all parameterizations of the bulk optical properties, the greatest density of the largest errors occurs in the far-infrared regions, where the infrared emission will have the lowest values and so will not contribute significantly to the OLR flux from the cloud, and in any case for qi < 10−6 kg kg−1, the cloud optical depth will be very small at those wavelengths and qi values.
The normalized probability distribution functions (PDFs) of ε(λ), calculated for each of the three bulk optical properties, are presented in Figs. 7a–c. The mean of the PDF of ε(λ) in Kext shown in Fig. 6a is −2.96%, but the standard deviation is ±14.7%. However, as previously stated, the experimental uncertainty in Kext is about ±50%, and since Kext is the first-order term for radiative transfer, then the σ value in the Kext relative error is well within the experimental uncertainty. Moreover, the parameterization is within ±7.5% for about 92% of all the cases considered.
The normalized PDF of ε(λ) found for ω0 is shown in Fig. 7b, and, as can be seen from the figure, the error is broader than Fig. 7a (not surprisingly, since the more difficult ice absorption regions are included in the parameterization). The mean of the ω0 PDF is 14.02% with a σ value of 14.7%, with 74% of the cases being within ±12.5%. The bias in the relative error calculated for ω0 is positive, with a value of 0.031. A positive bias, with respect to the true value of ω0, could nudge the cloud to greater absorption (because of the predicted value of ω0 being lower than the true value) and, therefore, warming, which could affect the temperature structure of the troposphere. However, this effect will be smaller than the impact of the first-order term Kext.
The distribution of the normalized PDF found for g is shown in Fig. 7c, and the figure shows that the parameterized g is within ±2.5% of the true value of g for about 80% of all values. Moreover, there is no significant bias and the mean of the PDF is −2.133% with a σ value of ±16.8%. However, this spread in g is primarily due to the longer wavelengths in the far-infrared, as shown in Fig. 6c, which in the PDF occur very infrequently. The impact of the parameterizations on the global model is discussed in the next section.
5. The impact of the parameterizations on GA5: A series of experiments
In this section, the impacts of the parameterizations on the GA5 configuration of the global model are presented in the form of a series of experiments under various assumptions. As previously stated in the introduction, each of the experiments is run over a 20-yr period, starting in 1989 and finishing in 2008. Each of the experiments is compared against the latest CERES annual 20-yr averaged shortwave and longwave mean fluxes at TOA, as compiled by Stephens et al. (2012), which uses the CERES reanalysis performed by Loeb et al. (2009). Each of the experiments is also compared against the ERA-Interim analysis of atmospheric temperature observations covering the same 20-yr period (Dee et al. 2011).
The control run assumes the following ice microphysics and ice optical parameterization. The mass–dimensional relationship is assumed to be the one from Brown and Francis (1995), and the PSDs are derived from the Houze et al. (1979) parameterization. The ice optical parameterization used in the control run is taken from E07. Each of the experiments is described in the next few paragraphs:
Experiment 1: Consistent PSD and inconsistent mass–dimensional relationship. The ice microphysics assumed in this experiment is taken from Furtado et al. (2014), and the F07 PSD parameterization is used for the ice crystal size spectrum. The ice optical parameterization is taken from Tables 1 and 2. In the radiation scheme, the F07 PSDs are generated using the mass–dimensional relationship predicted by the ensemble model (i.e., mass = 0.04D2). The bulk mass extinction and scattering coefficients are determined by applying different weights to each member of the ensemble model for each of the size bins in the PSD. The weighted-averaged mass extinction and mass scattering coefficients 〈Kext/sca〉 are therefore given by the following equation:
where Dmin and Dmax are the minimum and maximum ice crystal dimensions in the PSD, n(D) is the F07 PSDs, and wj is the weight applied to the jth ensemble member, where by definition Σwj = 1 for each of the size bins in the PSD. In Eq. (9) the subscripts ext/sca have been removed from kj(D) for reasons of clarity. Obviously, Kext/sca are the mass extinction and mass scattering cross sections predicted for each ensemble member of maximum dimension D. The weights in this experiment are assumed to be 0.90, 0.07, and 0.03, applied to the first three ensemble model members, respectively. These weights are chosen because they are the same as those used by Baran et al. (2014a) in evaluating the ensemble model against solar and infrared measurements obtained above optically thin cirrus. However, it is as yet unclear if these weights are applicable over a broader range of optical depth than was sampled by Baran et al. (2014a).
Experiment 2: Consistent PSD and mass–dimensional relationship. In this experiment, the same ice microphysics is assumed as in experiment 1, except that in the radiation scheme, the F07 PSDs are generated assuming the Cotton et al. (2013) mass–dimensional relationship (i.e., mass = 0.0257D2), and these PSDs are used to predict the bulk ice optical properties. The weights applied in this experiment to obtain the weighted-average mass extinction and scattering coefficients are the same as those assumed in experiment 1.
Experiment 3: Fully consistent but ensemble model reweighted toward more ice aggregated members. This experiment assumes, in the radiation scheme, the same mass–dimensional relationship as used in experiment 2, and the ice microphysics is kept the same as in experiment 1, except that in this experiment the ensemble model mass extinction is weighted toward the more aggregated members of the ensemble. The weights assumed here are 0.1, 0.2, 0.1, 0.4, 0.1, and 0.1, respectively.
Experiment 4: Fully consistent but ensemble model reweighted toward less aggregated members. In this experiment, the control run has been changed. Here the control ice microphysics is the same as in experiment 1, but the radiation scheme, which is the E07 parameterization, is kept the same as the previous control. In the experimental run, the F07 PSDs used in the radiation scheme have been generated assuming the same mass–dimensional relationship as used in experiment 2, but the ensemble model is weighted so that the mass extinction coefficients are equivalent to the mass extinction coefficients of experiment 1. To obtain this equivalence, the ensemble model was weighted using the following set of weights: 0.5, 0.2, 0.3, 0, 0, and 0. Here, the PSDs and mass–dimensional relationships assumed in the cloud and radiation schemes of GA5 are completely consistent. This experiment is hereafter referred to as the “consistent radiation treatment,” as this experiment is consistent with the solar and radiometric observations of Baran et al. (2014a) but using the Cotton et al. (2013) derived mass–dimensional relationship to generate the PSDs. The results obtained from each of the experiments are now discussed in sequential order.
a. Experiment 1
As a reminder, in this experiment, the new ice microphysics is utilized, and the ensemble model mass–dimension relationship is used to generate the F07 PSDs. The shortwave and longwave cirrus radiative properties are calculated using the new ice optical parameterization given in Tables 1 and 2, respectively. The shortwave and longwave flux differences at TOA between the simulations and the CERES observations are shown in Figs. 8a–d. As shown in Fig. 8b, experiment 1 leads to a slight improvement in the area-averaged SW root-mean-square error (rmse) of 0.13 W m−2. However, there are discernible divergences in the spatial distribution of the differences between the control and experiment 1. For instance, in the tropics, experiment 1 is not as bright as the control, especially off the coast around Brazil and off the West African coast near the equator, as well as centered on the equator over the Indian Ocean. Note also, experiment 1 over Antarctica is significantly less bright than the control. There is also a distinct decrease in the shortwave reflection, relative to the control, around the west Atlantic Ocean to near-neutral values and around the UK. The OLR differences are shown in Figs. 8c and 8d.
In the longwave, the new parameterization improves the area-averaged rmse by 0.21 W m−2. Moreover, off the coast of Antarctica, the new parameterization decreases the OLR differences to generally neutral values, although in places, such as the high midlatitudes over the Atlantic Ocean and over Canada, experiment 1 decreases the cloud transmission relative to the control. In general, over the tropics, experiment 1 increases the OLR relative to the control by about 5–10 W m−2; therefore, the ice cloud in experiment 1 is generally more transmitting in that region relative to the control. With respect to the control, the new parameterization does not have any undesirable impacts on the radiative properties of the model. However, the new parameterization has its biggest impact on the model in terms of the temperature structure of the troposphere, as shown in Fig. 9b. The figure shows that, in the tropical troposphere, experiment 1 removes the cold bias that is evident in the control shown in Fig. 9a, probably because the ice cloud in the control, over the western Pacific, is more transmitting relative to the cloud in experiment 1. However, experiment 1 does produce a slight warm bias of about 0.5 K in the southern midlatitudes. In general, in the climatically important tropical region, experiment 1 and the new ice optical parameterization clearly have a positive impact. The impact of changing the mass–diameter relationship in the radiation scheme is discussed in the next experiment.
b. Experiment 2
In this experiment, the ice microphysics is the same as experiment 1, but the F07 PSDs in the radiation scheme are generated using the Cotton et al. (2013) mass–dimensional relationship, and these PSDs are used to predict the ice bulk optical properties. The impact of experiment 2 on the SW and LW flux differences at TOA is shown in Figs. 10a and 10b, respectively. The impact of experiment 2 on the temperature structure of the troposphere is shown in Fig. 11. Clearly, the figures show that, on broadening the PSD in the radiation scheme while keeping the size–shape distribution the same as in experiment 1, the impact on the radiative fluxes and temperature structure of the troposphere is clearly considerable. In this experiment, relative to experiment 1, the area-weighted SW and LW rmse have increased by 0.68 and 1.29 W m−2, respectively. Both Figs. 10a and 10b show that, on broadening the PSD, there is considerably more SW flux reflected back to space and considerably less OLR relative to experiment 1. In general, it can be concluded from experiment 2 that the ice cloud is generally less transmitting than the ice cloud predicted in experiment 1. The impact of less transmitting cloud on the temperature structure of the troposphere is demonstrated in Fig. 11. The figure shows that the broadening of the PSD, which allows the occurrence of larger ice crystals relative to those in experiment 1, results in more absorption and therefore more warming of the troposphere at high altitudes in the tropics by about 0.5 K and by about 1 K in the southern midlatitudes. There is also some warming, of about 0.5 K, at altitudes greater than 600 hPa at northern midlatitudes.
Clearly, this experiment demonstrates the need to constrain the general shape of the PSD and determine more accurate mass–dimensional relationships. The impact of changing the ensemble shape distribution while preserving the shape of the PSD to be the same as that used in experiment 2 is discussed in the next experiment.
c. Experiment 3
For this experiment, the same ice microphysics and PSDs are used as in the previous experiment, but the ensemble model is reweighted toward the more ice aggregated members to increase ice cloud transmission. The results of experiment 3 for the SW and LW and its impact on the temperature structure of the troposphere are shown in Fig. 12a, Fig. 12b, and Fig. 13, respectively. Relative to experiment 2, increasing the ice cloud transmission by reweighting the ensemble model to more aggregated members results in less bright cloud and greater transmission in the OLR, as shown in Figs. 12a and 12b, respectively. Some of the SW impacts are desirable, especially over the Atlantic Ocean, where differences between the model and CERES are close to neutral values. This is also true in certain parts of the Southern Ocean but is by no means general. Interestingly, over Antarctica the impact of experiment 3 is similar to that of experiment 1, and both lead to an improvement over that region relative to the control. However, in general, experiment 3 clearly increases the area-weighted rmse by 1.7 W m−2 relative to experiment 1, which is undesirable.
The impact of increasing the cloud transmission, relative to experiment 1, is to reintroduce the cold bias seen in the control in the tropical troposphere, as shown in Fig. 13. It is also clear that experiment 3 leads to an increase in the cold bias in the lower atmosphere at high northern midlatitudes relative to the control and experiment 1. This experiment demonstrates that changing the values of the weights applied to habit mixture models of cirrus will also have an important impact on the radiative fluxes at TOA and on the temperature structure of the atmosphere. Therefore, it is of necessity to constrain not only the shape of the PSD but also the size–shape distribution across the PSD. Experiments 2 and 3 demonstrate that the shape of the PSD and the weighting of habit mixture models are of equal importance.
d. Experiment 4
In this last experiment, the F07 PSDs are the same as those used in experiment 2, but the ensemble model is reweighted so that the ice mass extinction and scattering coefficients are the same as those used in experiment 1. In the case of the control run, the ice microphysics is the same as was used in experiment 1, but the ice cloud radiation scheme is kept the same as in the previous control run. Figure 14 shows the reflected shortwave (Figs. 14a,b) and OLR (Figs. 14c,d) differences at TOA between E07 and the new ice optical parameterizations and the CERES observations. The figure shows that the consistent radiation treatment produces less bright cloud and more longwave cloud transmission relative to the E07 parameterization, which is an expected result from Figs. 5a and 5b. However, the bias toward more longwave cloud transmission results in more areas of near-zero differences between model and observation relative to the E07 parameterization. This is especially true in the Southern Ocean, and areas in the midlatitudes, as shown in Fig. 14d. The results obtained from this consistent radiation treatment are very similar to the results from the inconsistent experiment 1 (Figs. 8b,d with SW rmse = 8.61 W m−2 and LW rmse = 6.97 W m−2).
To understand these differences between the two radiation schemes and the CERES observations, Figs. 15a and 15b show the SW and LW flux differences between the consistent radiation and the E07 parameterization. Figure 15a shows that most of the SW differences between the two radiation schemes occur in the tropics, where the consistent radiation treatment predicts cloud that is generally less bright than the E07 parameterization by about 10 W m−2, and in the tropical warm pool region this difference can be up to about 15 W m−2. The outgoing LW differences between the two schemes are shown in Fig. 15b, and from this figure, the consistent radiation treatment allows more cloud LW transmission than the E07 parameterization. Outside of the tropics, the differences are up to about 5 W m−2, while in the tropics the differences can be up to about 15 W m−2. The higher LW differences are generally around the tropical warm pool.
Given these SW and LW differences between the two radiation schemes and the fact that the consistent radiation is a coupling between the microphysics and bulk optical properties, it would be interesting to compare differences in the SW and LW radiative effect between the two schemes and the CERES observations. The SW and LW radiative effect differences between the model and CERES are shown in Figs. 16a–d, assuming the E07 parameterization and the consistent radiation scheme, respectively. Figure 16a demonstrates that, given improved microphysics, the SW radiative effect that the E07 parameterization predicts is sometimes overcooling relative to the CERES product. In the tropics and midlatitudes, these differences can be up to about −15 W m−2, and in some tropical locations the radiative effect differences can reach up to about −30 W m−2 off the coast of South America and along the equator above the Indian Ocean. The SW radiative effect predicted by the consistent scheme is shown in Fig. 16b, and this essentially shows the opposite to the E07 parameterization. The consistent scheme generally decreases the cooling predicted by the E07 parameterization, although relative to CERES, this can appear as too much warming by up to about 30 W m−2 in some locations in the tropics and especially around the tropical warm pool. However, in the midlatitude regions (with the notable exception of the tropical warm pool), the consistent radiation scheme appears to decrease the cooling predicted by the E07 parameterization, which appears desirable relative to the CERES product.
The LW radiative effect predicted by the E07 parameterization and the consistent radiation scheme is shown in Figs. 16c and 16d, respectively. The two schemes, relative to the CERES product, generally appear more similar to each other. However, in the tropics and midlatitudes, the consistent radiation scheme tends to offset the E07 LW radiative effect more in the direction of the CERES observations. This similarity in the longwave radiative effect between the two schemes could be because the model predicts the cloud to be at certain altitudes, and it is the radiating temperature of these clouds at those altitudes that dominates the CERES measurements, rather than differences in the radiation parameterizations. The differences in the SW and LW radiative effects between the consistent radiation scheme and the E07 parameterization are shown in Figs. 17a and 17b, respectively. Figure 17a shows that the consistent radiation scheme generally predicts less SW cooling than the E07 parameterization by about 5 W m−2 outside of the tropics; in the tropics this can reach up to about 15 W m−2, while over the tropical warm pool, the difference can reach up to about 20 W m−2. Figure 17b shows that the differences in the LW radiative effect between the two radiation schemes occurs either in the tropics, generally around the tropical warm pool, or at high latitudes in both hemispheres, and the signs are opposite. Not surprisingly, given previous results, over the tropical warm pool the consistent radiation predicts less LW radiative effect than the E07 parameterization by about −10 W m−2, while at the high latitudes in both hemispheres, the consistent radiation scheme predicts more of an LW radiative effect than the E07 parameterization, by up to about 10 W m−2.
The results presented in this section have shown that, given any mass–dimensional relationship, the equivalent fluxes, or for that matter, radiances can be predicted by reweighting the ensemble model, given observational radiometric constraints on the distribution of weights or in situ estimates of the size–shape PSD.
The impacts of the E07 and consistent radiation parameterizations on the temperature structure of the troposphere are shown in Figs. 18a and 18b, respectively. Figure 18b shows that the consistent radiation treatment is generally similar to experiment 1 (Fig. 9b), most notably, with a further slight reduction of the cold bias in the tropical troposphere relative to Fig. 18a. This same cold bias is also significantly reduced relative to the original control by the E07 parameterization when convolved with the more representative cloud microphysics derived from Furtado et al. (2014) (Fig. 18a). There is a further slight warming of southern midlatitudes in the troposphere produced by the consistent radiation treatment relative to the E07 parameterization. Differences in the zonal mean temperature structure of the atmosphere between the consistent radiation scheme and the E07 parameterization are shown in Fig. 19, and this figure shows that the consistent radiation scheme is generally more warming than the E07 parameterization by up to about 1 K.
This experiment has shown that the consistent radiation treatment predicts cloud that is generally less bright and increases the OLR relative to the E07 parameterization, and generally warms the zonal mean temperature structure of the atmosphere by up to about 1 K. These effects have the consequence that the consistent radiation scheme will predict less of an SW radiative effect than the E07 parameterization. However, differences in the LW radiative effect between the consistent radiation scheme and the E07 parameterization mostly occur around the tropical warm pool, where they are negative, and at high latitudes in both hemispheres, where the differences are positive. Therefore, given the same microphysics scheme in GA5, the consistent radiation scheme predicts less of an LW radiative effect around the tropical warm pool but more of an effect at high latitudes in both hemispheres. A comparison between the CERES global means predicted by the consistent radiation treatment and the E07 parameterization is shown in Table 3. The Table shows that the E07 parameterization has increased the outgoing shortwave and decreased the outgoing longwave fluxes at TOA relative to the consistent radiation treatment, in accord with expectation. This is most notable in the shortwave, where the E07 parameterization now exceeds the observational uncertainty in the outgoing shortwave flux at TOA. It is expected that given further improvements to the global model, which may result in more ice mass being predicted, the E07 parameterization will result in greater divergences between model and observation. However, Table 3 also shows that the consistent radiation treatment in the absorbed shortwave and SW CRE are both outside the observational uncertainty. It is expected that, with an improved climate model prediction of ice mass, these two SW radiative components should agree within the observational uncertainty. A further possibility as to why the consistent radiation treatment underpredicts some of the SW components shown in Table 3 could also be due to the asymmetry parameter predicted by the ensemble model not being sufficiently low. However, for now, it is more likely that improving the representation of ice mass in the Met Office GA configuration will improve the consistent radiation prediction of those two SW components. Although the E07 parameterization is not detrimental to the performance of GA5 with improved microphysics, with further improvements to the model, the E07 parameterization can only get worse when compared against observation. This is because the cancellation of error that was discussed in the introduction will be removed through further model improvements. Consistent physics within predictive models must always be preferred over inconsistent treatments. Unfortunately, the inconsistent treatments are currently prevalent in climate models, which may give reasonable predictions, but for the wrong physical reasons.
6. Conclusions and discussion
A new flexible shortwave and longwave parameterization of the ensemble model bulk optical properties has been presented and applied to the Met Office Unified Model Global Atmosphere 5.0 configuration (GA5). The bulk optical properties have been computed between the wavelengths of 0.2 and 120 μm at sufficient resolution in wavelength space such that the rapid variation of the imaginary refractive index of ice in the near-infrared and longwave parts of the spectrum is captured. It has been shown that the parameterization can reproduce the database of 2.99 × 106 bulk mass extinction coefficients to within ±7.5% for over 90% of the database. The largest errors occur at wavelengths in the far-infrared. At these wavelengths, the Planck function will be small and so will not bias the values of mass extinction. Currently, the mass extinction of cirrus is not known experimentally to within ±50%. The mass extinction is the first-order term that determines the reflection and transmission properties of cloud. The error in the mass extinction parameterization is well within the current experimental uncertainty. The single-scattering albedo and asymmetry parameter parameterizations were found to be within ±12.5% and ±2.5% for 74% and 80% of the database, respectively. The largest errors in the single-scattering albedo and asymmetry parameter tended to occur at the lowest values of IWC and occurred in the far-infrared regions.
The new parameterization is more flexible than previous parameterizations used in the Met Office series of GA models, since the weighting of the ensemble members can be changed, given some mass–dimensional relationship to generate the parameterized PSD, and this preserves physical consistency between the cloud physics and radiation schemes of the climate model. This flexibility was demonstrated in experiments 3 and 4. Experiments 2, 3, and 4 demonstrate that the weighting of habit mixture models and assumed shapes of the PSDs or mass–D relationships are of equal importance.
The paper has demonstrated that the overarching philosophy in demanding consistent cloud physics between the cloud and radiation schemes does not have detrimental effects on the overall performance of GA5. Rather, if this overarching philosophy is not followed, as demonstrated by the control used in experiment 1 and the application of the E07 parameterization used as the control in experiment 4, then inconsistent microphysics may lead to desirable results, but for perhaps the wrong physical reasons. Insisting on consistency between the cloud and radiation schemes means that all other parameterizations that are predicting the cloud macrophysical and microphysical properties must realistically represent the most pertinent processes relevant to those parameterizations.
The experiments presented in this paper demonstrate the importance of observationally constraining the shape of the PSD and which general mass–dimensional relationship to apply to climate models. Furthermore, just as important as this is to determine experimentally which size–shape distribution (in terms of a weighted habit mixture model of cirrus) best matches global SW and LW radiometric observations.
Relying on error cancellation, as discussed at the beginning of this paper and in the discussion of experiment 4, is undesirable, as when the distributions and amounts of IWC are finally addressed in climate models, those parameterizations that do rely on error cancellation can only get worse (assuming IWC is increased to the observed amounts and at the correct altitude). This, to some extent, has already been identified in the results presented in Table 3. A further overarching philosophy of the approach in this paper is to remove error cancellation by directly relating global model prognostic variables to ice optical properties. In this way, global model prognostic variables are directly tested against radiative measurements, rather than placing a reliance on diagnosed variables, which may lead to error cancellation; this will then preclude identification of systematic climate model biases.
We thank Professor Greg McFarquhar and two anonymous reviewers for their helpful comments.