A reduced-form approach to characterizing sulfate aerosol effects on climate in integrated assessment models. Final report

The objective of this study was to devise a methodology for estimating the spatial patterns of future climate change accounting for the effects of both greenhouse gases and sulfate aerosols under a wide range of emissions scenarios, using the results of General Circulation Models.


OBJECTIVE:
To devise a methodology for estimating the spatial patterns of future climate change accounting for the effects of both greenhouse gases and sulfate aerosols under a wide range of emissions scenarios, using the results of General Circulation Models.

BACKGROUND:
The realization that anthropogenic aerosol effects on climate may be as important as those of C02 and other greenhouse gases has radically changed our ideas regarding how climate may change in the future. Historically, the potential importance of aerosols has been known for decades (see, e.g., Mitchell, 1975). What has changed in recent years has been our ability to quantify the direct (clear-sky) and indirect (cloud albedo) effects. Many important papers have been published on: the physical principles (e.g., Twomey et al., 1984;Charlson et al., 1991); feedbacks that may enhance the effects of natural precursor emissions (Charlson et al., 1987); global-scale climate consequences of anthropogenic SO2 emissions (Schwartz, 1988;Wigley, 1989Wigley, ,1991; patterns and magnitudes of radiative forcing effects (Charlson et al., 1991;Kiehl and Briegleb, 1992;Jones et al., 1994); detailed climate response (Taylor and Mitchell et al.;1995;Roeckner et al., 1995;Erickson et al., 1995); and consequences for detection of human-induced climate change (Santer et al., 1995b(Santer et al., ,1996aKarl et al., 1995). Reviews have been published by, among others, IPCC (Shine et al., 1995(Shine et al., ,1996Jonas et al., 1995Jonas et al., , 1996, Charlson et al., (1992), Penner et al. (1994), Charlson and Wigley (1994)-see also Wigley (1994Wigley ( ,1995. The importance of sulfate aerosols is two-fold: they are crucial in explaining past climate change; and in predicting future change. For the former, incorporation of aerosol effects has allowed IPCC to go from a 1990 position where it was noted that we could not identify the spatial signature of anthropogenic climate influences (Wigley and Barnett, 1990) to one where these influences have been identified with a high degree of statistical confidence (Santer et al., 199610). For prediction, the most recent coupled ocean-atmosphere GCM studies at the Hadley Centre, UK, and the Max-Planck Institute, Germany, show that the regional details of future climate differ radically from earlier C02-only projections when the influence of S02-derived aerosols is included in the calculations (Kattenberg et al., 1996).
The influence of aerosols on future climate poses a serious problem for Integrated Assessment modeling. For the major part of the forcing due to greenhouse gas concentration changes, the patterns of future climate change do not depend on the geographical distribution of emissions. This is largely because gases like C02, C-, N20 and many halocarbons are long-lived and, hence, are well mixed in the atmosphere. For aerosols and their precursors, however, their short lifetimes (a few days) mean that they have spatially heterogeneous forcing effects that depend on the locations of emissions. This produces spatially heterogeneous climate response patterns that also depend on the emissions patterns. For the climate response to SO2 emissions, therefore, we need to know not only when, but also where the emissions will occur. Thus, for each future emissions pattern possibility, a separate climate change simulation must be carried out-unless some simpler method can be devised.
Fortunately, a relatively simple method does exist, described below. Its basis is the scaling method for developing regional climate change scenarios devised by Santer et al. (1990)-see also Hulme et al. (1994Hulme et al. ( ,1995. This method was devised to allow GCM results to be extrapolated easily to a wide range of greenhouse-gas emissions scenarios, and to allow uncertainties to be assessed. The method is to decouple the global-mean and regionalization components of future climate change using where AT(t) is the global-mean temperature change and Ak($ is the normalized pattern of change (i.e., per degree global warming) for climate variable "X'. AT(t) can be obtained from an upwelling-diffusion, energy balance model-still the main tool used by IPCC for such projections (Kattenberg et al., 1996). AT(t) simulations may be made quickly for any given emissions scenario, while uncertainties can be quantified by changing the UD-EBM parameters (most importantly, the climate sensitivity ATh).
k ( x ) information may be obtained from either: equilibrium, mixed-layer ocean (ML0)-atmosphere GCM experiments; or coupled ocean-atmosphere GCM experiments (provided a sufficient change has occurred to define the signal adequately relative to the model's internal-variability noise). Additional uncertainty information may be obtained by examining inter-model differences in &(IC), and uncertainties may be reduced by averaging results from different models.
The validity of equ.
(1) depends on whether the signal pattern, as it evolves under constant forcing, remains constant in time. As far as we can judge, this is a correct assumption (see, e.g., Santer et al., 1995a). Use of MLO-AGCM results is valid provided their signal patterns agree with those of 0-AGCMs. This is also a justifiable assumption except in regions of oceanic deep-water production (where no one lives).
In the following, the UD-EBM employed is MAGICC (Model for the Assessment of Greenhouse-gas Induced Climate Change), which is the same model employed by IPCC. Further details of MAGICC are given in Raper (1987,1992), Raper et al. (1996) and Hulme et al. (1994). The regionalization software used is SCENGEN (Hulme et al., 1995), currently being updated to employ the results of this project.

RESULTS:
The results described below are divided into three sections: modifications to MAGICC; algorithm development for SCENGEN; and application to MERGE.

I
This work has been supported jointly by the present project and by funds from the National Science Foundation.
When it was decided to use MAGICC for IPCC calculations of future global-mean temperature and sea-level projections, the "core" UD-EBM had to be ungraded to incorporate the latest results from 0-AGCMs and validated against such models. The main modifications were to account for the fact that the climate sensitivity over ocean differs from that over land, and the fact that the intensity of the thermohaline circulation (or upwelling rate, w, in the UD-EBM) apparently slows down as the globe warms. To include differential land/ocean sensitivity effects we required the input parameters to be the global sensitivity and the land/ocean sensitivity ratio. Since the UD-EBM requires information on the feedback parameter values over land and ocean in each hemisphere, some additional mathematical analysis was required to relate the two parameter specification methods. The variable-w aspect had been built into the model initially-we simply added a feedback loop so that "w" was a function of global-mean ocean temperature. Further details are given in Raper et al. (1996) and Kattenberg et al. (1996). In the proposal for this project, it was thought that the UD-EBM may have to be modified to account for a possible dependence of the climate sensitivity on the type of forcing. The original work of Taylor and  seemed to indicate that the sensitivity to aerosol forcing (i.e., the equilibrium global-mean warming for unit radiative forcing) may differ from that for C02. It was subsequently found  that this result was an artifact of the GCM not having reached equilibrium. When longer runs were completed, there was no significant difference between the C02 sensitivity and the aerosol sensitivity.

Algorithm Development for SCENGEN
The method that is currently being implemented is as follows. The basic assumption is that the patterns of climate-change response from different forcings can be linearly superimposed to give the combined response. This is a standard assumption that has always been made, albeit tacitly, for greenhouse-gas forcings (and emissions). It is the basis for the "equivalent-C02" concept, and is justifiable in this case because the most important greenhouse gases are well-mixed in the atmosphere. For combined C02aerosol forcing, however, the assumption may not be correct because of the strongly regional nature of aerosol forcing (in turn, a consequence of the short lifetimes of aerosols and their precursors). We have therefore tested this assumption extensively for temperature  using the results of Taylor and . It works extremely well.
We have, as part of this project, carried out further tests for precipitation. These results were inconclusive, however, because the signal-to-noise ratios of precipitation responses in the Taylor and Penner results are too low to be able to detect significant differences between the response patterns for combined C02-aerosol forcing, and the sum of the individual C02-and aerosol-forcing response patterns. Currently, there are no suitable model experiments available to test the superposition hypothesis adequately for variables other than temperature.
The superposition method is expressed by the equation: where AkC(x) is the normalized response pattern for variable X for forcing by gases that give a CO2-like response, and d a , i ( x ) is the normalized response pattern for the i-th SO2 emissions pattern. ATc(t) and ATa,i(t) are the corresponding global-mean temperature changes (from MAGICC). Equation (2) is a straightforward extension of the Santer et al. (1990) method (equ. (1)).
In the above spatial breakdown, we separate the aerosol and non-aerosol components of global-mean temperature change, since they have different spatial patterns. Ideally, one should separate all factors that have different spatial patterns.
We make an assumption here that the signal can be divided into just two basic types of response pattern: a CO2-like pattern and a set of aerosol-like patterns (cf. Penner et al., 1996). We lump all greenhouse gases together as "C02-like" and both direct and indirect aerosol effects together as aerosol-like (strictly, sulfate-aerosol-like).
To apply this equation, we require that any future SO2 emissions data are given separately for each of I regions (i = 1, I): currently we use only three regions, a North Atlantic region (Europe/North America), an Asian region (India, China, Japan), and the rest of the world. Climate response pattern data are not yet available for these canonical cases, but should become available within the next year or so. The ATc and AT, scaling factors come directly from MAGICC by running the model with the I + 1 forcing cases separately.
There are alternatives to equ. (2), but we judge these to be less suitable for the MAGICC/SCENGEN case. For example, as suggested in the original proposal, it would be possible to use the component radiative forcings as scaling factors for the canonical patterns. The two methods reduce to the same at equilibrium, but the temperature weighting method is clearly better in the time-dependent case because it accounts for the non-linear relationship between forcing and (lagged) response.
In order to use equ. (Z), we need to specify the forcing time series with which to drive MAGICC in order to determine the temperature weights, ATC(t) and ATa,i(t). This is straightforward for ATJt), but less so for the aerosol weighting factors. To explain why, we need to give more detail on how MAGICC specifies aerosol radiative forcing.
MAGICC does this using simple relationships between forcing (Q) and SO2 emissions (E) given by Wigley (1991); see below. The direct (clear-sky) forcing is a linear function of emissions level, while the indirect (cloud-albedo) forcing is logarithmic in emissions. Because the latter forcing is non-linear in emissions, deriving its regional decomposition is a non-trivial problem: for any given Q-E relationship, forcing for the sum of the regional emissions will not be equal to the sum of the individual forcings (see Figs. 1 and 2: the full explanation for these Figures is given below). Hence, some form of scaling of the regional indirect forcings is required to ensure consistency with the global forcing as determined by the global emissions level. It is scaling to the correct global sum that is necessary, rather than modifying the global  (Leggett et al., 1992) as a baseline, extended to 2200 by assuming that these emissions decline to a value 10Tg S/yr below the current level (i.e., to 65Tg S/yr) by 2200. The regional breakdown is highly idealized, chosen for illustrative purposes only. Region 1 (North Atlantic) is assumed to have a continuous decline in emissions-as has occurred in recent decades-to a total value close to zero in 2050. Region 2 (Asia) is assumed to be the main source region for future SO2 emissions, while Region 3 (Rest of the World) is assumed to have emissions that constitute 10% of the total fof Regions 2 and 3. This total in turn is derived by subtracting Region 1 emissions from the extended IS92a global emissions. YEAR Figure 2 Global sulfate aerosol radiative forcing breakdown (direct plus indirect effects) for the SO2 emissions scenarios given in Fig. 1. The sum of the independentlycalculated individual regional forcing results can be seen to be significantly smaller than the forcing obtained using total global emissions.
sum, because the primary determinant of the magnitude of indirect forcing is its 1990 global value (see Raper et al., 1996;Kattenberg et al., 1996).
To explain how this scaling is done we introduce the following terminology.. . .... ( 5) where Enat is a natural SO2 emissions level taken as 42Tg S/yr and AQ1(1990) = -0.8W/m2 is the value used by IPCC. The 1990 AQ values are highly uncertain, especially the indirect value, so are left as free parameters in the code. In Wigley and Raper (1992), and in IPCC, a slightly different formulation is used to account for changes in SO2 emissions heights (and, hence, lifetime), but this secondary effect does not seem warranted here at present. We assume initially that the same formulae may be applied at the regional level. This is perfectly satisfactory for the direct forcing, but requires an adjustment for the indirect forcing.
Up to 1990, MAGICC is forced by the sum of the forcings due to the various greenhouse gases, a biomass-derived aerosol term, and sulfate aerosols according to equs. (4) and (5) Raper et al., 1996, for further details). After 1990, to obtain the global-mean temperature changes associated with each separate SO2 emissions region (ATa,i(t)) MAGICC must be forced by the non-sulfate aerosol forcing term and by each regional SO2 emissions case separately (i.e., four cases). In each case, full forcing must be used to 1990 (i.e., full global SO2 emissions ) and either zero or single-region SO2 emissions subsequently. Thus, the required S@ input for MAGICC is the global SO;! emissions history to 1990, and the regional changes in SO2 emissions relative to 1990 after this date. The regional sulfate aerosol forcing is therefore given by .. . . (7) Before doings so, it is useful to examine the magnitude of the nonlinearity effect.
From equ. (5) it can be shown that the ratio of the sum of the independently-calculated regional indirect forcings to the indirect forcing based on global emissions is given by As an example, suppose C6Ei = 75Tg S/yr and consider two different regional breakdowns. First, R would be 1.17 if 6Ei =25Tg S/yr for i = 1,3; while R is 0.76 if 6E1= -25Tg S/yr, 6E2 = 100Tg S/yr and 6E3 = OTg S/yr. Note that the ratio may be greater or less than one. The example shown in Fig. 2 is qualitatively similar to the second of the above simpler examples (i.e., R e 1.0). Figure 2 shows that the forcing error that would arise by summing the unmodified results for individual regions is likely to be quite small (less than O.lW/mZ in this case, which is much less than the total forcing). Thus, ignoring this nonlinearity problem would lead to only small errors in future temperature projections. More importantly, errors that might arise from using an incorrect method to apportion the regional forcings would be very small indeed. We now return to devising an appropriate method.
The sum of equ. (7) over i obviously will not give the correct total temperature response since each individual case repeats the pre-1990 changes. This is not an issue because SCENGEN uses only the incremental temperature changes from 1990 (or some earlier reference period, such as 1961-90, for which a baseline climatology is available). What is of concern here is that the sum over the regions of the post-1990 incremental forcings (and, hence, temperature changes) will not equal the single globally-forced value. In other words, while . . . .( 11) even though it produces the correct forcing sum, will not work. It clearly fails if CGQ1,i = 0, which is quite possible, and it cannot give the correct limiting result if emissions were to return to pre-industrial levels.
Instead, we scale the total regional forcings from 1765 rather than from 1990. The sQ;,i(t) are then determined by To apply equ. (12), we need to know AQ1,i(1990); i.e., the regional indirect forcing breakdown in 1990. A reasonable approximation (to which, as noted above, the results are relatively insensitive) is to assume AQ1,i(1990) = PAQ1,g(1990 .... ( 13) where P = 0.5,0.5 and 0.0 for regions 1 (North Atlantic), 2 (Asia) and 3 (Rest of the World).
We have applied this method to detection studies in order to test the standard pattern correlation method of Santer et al. (1993Santer et al. ( ,1995a under controlled conditions . In these studies, we use equ. (2) to generate the time-depedent response to past C02-plus-aerosol forcing (the signal) to which we add members of sn ensemble of noise simulations (from the control-run of an 0-AGCM). There is scope to apply this noise-addition. method to future SCENGEN-type climate change scenarios to quantify uncertainties related to the natural variability of the climate system and/or to investigate some aspects of extreme event analysis.

Application to MERGE
Incorporation of aerosols into MAGICC and SCENGEN represents only a part of their inclusion into IA modeling in general, since MAGICC and SCENGEN do not comprise a full IA model. However, MAGICC and SCENGEN (at least their earlier versions) are included in the Battelle PNL integrated assessment models, and it is planned to update the PNL models when the present methodology has been fully implemented and appropriate GCM data are available.
In the meantime, it was considered worthwhile trying to account for aerosols in a much simple LA model, MERGE (Manne and Richels, 1992;Manne et al., 1995). The development below uses a differential equation formulation that differs from the MERGE finite difference description. The basic approach is as in MAGICC/SCENGEN, i.e., to consider the inclusion of sulfate aerosols at the global level first, and then to disaggregate the results to the nine regions used by MERGE. The global perspective follows MAGICC, except that the climate model is a one-box model as currently used by MERGE. Extension to a two-box model would be straightforward and would improve the realism of the climate model substantially. The spatial disaggregation of aerosol effects is based on the rather crude assumption that the regional transient temperature changes due to aerosol effects depend linearly on the regional emissions changes.
The global-mean model used by MERGE may be written as .... ( 14) where AQ is the anthropogenic radiative forcing change (ghg = greenhouse gases; dir = direct sulfate aerosols; ind = indirect (cloud albedo) sulfate aerosols; other = other anthropogenic forcing), C is a thermal inertia term, a defines the response time of the system, and all changes are measured relative to an initial preindustrial steady state (t = 0 is generally taken as 1765). In equilibrium, so that the equilibrium C02-doubling temperature ch&qge (AT2x; i.e., the climate sensitivity) is given by The product aC is therefore inversely proportional to the climate sensitivity.
In the MAGICC /SCENGEN case, we follow through four separate components of global-mean temperature change as pattern weighting factors (equ. (2) above). In the MERGE case, we only follow through two components, the CO2-like and aerosol-like parts. Thus, the global-mean climate model calculations need only to be carried out using equ. (14) and another equation that isolates the CO2-like component, given by This approximates the global-mean temperature change due to emissions with a C02like response pattern. Since the models (14) and (15) are linear, the aerosol-like component of global-mean temperature change is given by .... ( 16) The next task is to spatially disaggregate ATc&) and ATa(t) and then recombine them to obtain the regional temperature changes.
We denote the regions by subscript "i" (i = 1,9). Since the regions do not comprehensively cover the whole surface area of the globe, we may introduce a tenth region that represents the remaining area.
For the CO2-like component (ATc(t)) we can use the results from GCMs to determine how the regional changes relate to the global-mean change-this essentially follows the SCENGEN method. If ai is the temperature change in region i for a unit global-mean warming due to forcings with a CO2-like response pattern, then ATc,i(t) = aiATc(t) . . . . ( 17) MERGE would simply use a table of ai values from C02 simulations with GCMs to generate the various ATcti(t).
For the aerosol-like component (ATa(t)), disaggregation is a little more complicated. Further work is required to improve on the method developed here. We divide the aerosol cooling into a globally uniform term and a regionally-specific term related to SO2 emissions per unit areh. Emissions per unit area must be used since the same emissions spread over a wider area would lead to a lower mean aerosol concentration and, hence, a smaller cooling effect. Thus .... ( 18) where ATa,0 is the globally-uniform component, AEi is the total emissions change for region i, and Ai is the area of that region.
If equ. (18) is multiplied by Ai, then summation over i yields  ...( 19) where AE is the global SO;? emissions change (= CAEi) and A is the area of the global (= C u i ) . Equation (19) determines p, which can be substituted into equ. (18) to give . . . . (20) where is the fraction of the globe covered by region i. A more convenient form for equ. (20) obtains if we replace ATato by the fraction of ATa that this uniform temperature change term represents giving ATa,i(t) = ATa(t)[@ + (1 -@)(AEi/AE)/ai] The overall regional change in temperature is then the sum of equs. (17) and (22). Note that @ can be estimated quite simply from available GCM-based climate model results (such as those of Taylor and  by regressing the pattern of annualmean temperature response against the pattern of emissions. The constant term in such a regression analysis gives AT,,o directly, and the correlation coefficient provides an indicqtor of the validity of the regionalization approach.
Further work on implementing this method into MERGE is continuing in collaboration with Rich Richels and Alan Manne.