Portal de Periódicos da CAPES

The aerosol–Bénard cell effect on marine stratocumulus clouds and its contribution to glacial-interglacial cycles

2011; American Geophysical Union; Volume: 116; Issue: D10 Linguagem: Inglês

10.1029/2010jd014470

2156-2202

R. Z. Bar-Or, Hezi Gildor, Carynelisa Erlick,

Aeolian processes and effects

Journal of Geophysical Research: AtmospheresVolume 116, Issue D10 Climate and DynamicsFree Access The aerosol–Bénard cell effect on marine stratocumulus clouds and its contribution to glacial-interglacial cycles R. Z. Bar-Or, R. Z. Bar-Or Department of Environmental Sciences and Energy Research, Weizmann Institute of Science, Rehovot, IsraelSearch for more papers by this authorH. Gildor, H. Gildor Fredy and Nadine Herrmann Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem, IsraelSearch for more papers by this authorC. Erlick, C. Erlick caryn@dina.es.huji.ac.il Fredy and Nadine Herrmann Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem, IsraelSearch for more papers by this author R. Z. Bar-Or, R. Z. Bar-Or Department of Environmental Sciences and Energy Research, Weizmann Institute of Science, Rehovot, IsraelSearch for more papers by this authorH. Gildor, H. Gildor Fredy and Nadine Herrmann Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem, IsraelSearch for more papers by this authorC. Erlick, C. Erlick caryn@dina.es.huji.ac.il Fredy and Nadine Herrmann Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem, IsraelSearch for more papers by this author First published: 28 May 2011 https://doi.org/10.1029/2010JD014470Citations: 3AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinked InRedditWechat Abstract [1] Aerosol-cloud interactions, such as aerosol loading in convective clouds resulting in either precipitation suppression or cloud invigoration, in higher cloud tops, and in longer-lived clouds, are well known. Here we investigate a new aerosol-cloud interaction, the effect of aerosol loading on Bénard cells, on the stratocumulus cloud fraction, and ultimately on the climate over glacial-interglacial cycles, using a two-dimensional model running a million year continuous simulation. This radiative effect is observed only in marine boundary layer stratocumulus clouds that have a convective cellular structure. Recent research suggests that aerosols can switch the direction of convection in Bénard cells (from open cells to closed cells) by suppressing precipitation and therefore dramatically change the cloud fraction. The effect investigated in this work differs from previously known aerosol effects on convective clouds by its intensity and magnitude and has never been taken into account in past climate simulations. The results show that accounting for the aerosol–Bénard cell effect alone contributes a negative radiative forcing, affecting both the Northern Hemisphere mean annual surface temperature and ice volume. Adding the aerosol–Bénard cell effect to the direct radiative effect of dust and to the effect of dust on snow and ice albedo shows that the aerosol–Bénard cell effect plays a significant role in glacial-interglacial climate change, strengthening the earlier glacial cycles and creating a larger glacial-interglacial surface temperature amplitude while preserving the continental ice volume amplitude. Because of the model limitations, there are a number of uncertainties involved. However, the results serve to give a preliminary evaluation of the aerosol–Bénard cell effect at least qualitatively if not quantitatively. Key Points We investigated a new radiative forcing called the aerosol–Benard cell effect Alone the aerosol–Benard cell effect contributes a negative radiative forcing The total effect of dust is a larger glacial-interglacial temperature amplitude 1. Introduction [2] It has been shown that large fluxes of dust to the atmosphere, as observed during the ice ages, would have a strong climatic effect on glacial-interglacial cycles [Lambert et al., 2008; Winckler et al., 2008; Goudie, 2009]. The dust produces (1) a negative direct radiative forcing on surface due to the increase of the total aerosol optical depth (AOD), (2) a positive radiative forcing of dust on snow and ice surfaces due to the decrease in their albedo, and (3) other radiative forcings due to aerosol-cloud microphysical interactions. The role of the first two effects on glacial-interglacial cycles, specifically the potential for dust to affect glacial-interglacial cycles through a forcing of the global radiation budget, was investigated previously by Bar-Or et al. [2008]. In the work by Bar-Or et al. [2008], a cooling trend was found due to the direct radiative forcing by dust, and a significant impact on the speed of ice sheet retreat was found due to the effect of dust on snow albedo. Furthermore, it was shown that dust radiative forcing emphasizes the asymmetry in the glacial cycles. Regarding the third type of effect, past research has shown that high aerosol concentrations in convective clouds may lead to precipitation suppression and higher cloud tops [Ekman et al., 2007; Freud et al., 2008] and/or may affect cloud dynamical processes, leading to more invigorated and longer-lived clouds [Koren et al., 2005, 2010]. In the present work, we parameterize a microphysical aerosol effect on clouds related to Bénard convection cells, recently suggested to have a strong impact on climate dynamics [Albrecht, 1989; Kaufman et al., 2005; Stevens et al., 2005; Rosenfeld et al., 2006]. At present, climate models cannot resolve aerosol-cloud-climate feedbacks associated with changes in cellular convection in marine boundary layer stratocumulus, and here we demonstrate a first attempt to parameterize such a microphysical effect in a global model with long time scales. [3] Marine stratocumulus clouds in low wind shear environments commonly form Bénard convection cells, convection cells that appear spontaneously in a liquid layer when heat is applied from below [see Koschmieder, 1993]. Observations show two possible cellular structures, open Bénard cells and closed Bénard cells [Wood and Hartmann, 2006]. Open cells are formed where surface radiative heating is dominant, causing strong humid updrafts in the narrow cell borders and a weak downdraft of cool dry air in the cell center. Closed cells occur when radiative cooling at the top of the clouds takes over and causes a descent of cooled air in the cell borders, while weak humid updrafts in the center of the cell fill the missing air mass in the upper layers. [4] Rosenfeld et al. [2006] recently suggested a mechanism by which aerosols may cause open Bénard cells to close by suspending drizzling processes and enhancing radiative cooling at the cloud tops. According to Rosenfeld et al. [2006], high aerosol concentrations cause open cells to close, dramatically changing the oceanic cloud fraction (CF) from around 0.25 to almost unity [Wood and Hartmann, 2006; Wang and Feingold, 2009a, 2009b]. Schreier et al. [2010] found that in the case of ship tracks, the change in the direction of convection occurs rapidly over wide areas and may last from several hours to days, and therefore has potential to change the atmospheric radiation balance in a similarly dramatic fashion. Rosenfeld et al. [2006, p. 2503] stated that, "The proposed mechanism suggests that very small changes in the aerosols input to the marine boundary layer can have large impacts on the oceanic cloud cover and likely in turn on the global temperature, in ways that are not yet accounted for in the climate models." We label this previously unaccounted for aerosol-cloud microphysical radiative forcing effect "the aerosol–Bénard cell effect." Although the radiative effects of dust during glacial maxima, including "indirect" effects have been simulated [Takemura et al., 2009], the aerosol–Bénard cell effect has not previously been considered in climate models at all [Forster et al., 2007; Penner et al., 2001], and here we make a first attempt at modeling the effect using a reasonable sequence of parameterizations in order to evaluate its potential importance to glacial-interglacial cycles. [5] This paper is organized as follows. In section 2, we briefly review the model used and the new parameterization introduced to account for the aerosol–Bénard cell effect. In section 3.1, we investigate the aerosol–Bénard cell effect alone and show that it contributes a negative radiative forcing. In section 3.2, we investigate the aerosol–Bénard cell effect combined with the direct radiative effect of dust and with the effect of dust on surface albedo and show that the combination increases the glacial-interglacial surface temperature amplitude in the earlier cycles but has only a minor impact on the continental ice volume. A set of sensitivity tests to variations in a few critical parameters is presented in section 4, and we conclude in section 5. 2. Model and Methods [6] In this study, we use the LLN 2D NH model (Louvain la Neuve 2D Northern Hemisphere), a two-dimensional, latitude-altitude-time-dependent coupled zonally averaged climate–ice sheet model [Gallée et al., 1991, 1992; Loutre and Berger, 2000] built to simulate the Northern Hemisphere climate and ice sheet evolution on long time scales. LLN 2D NH simulates the Northern Hemisphere only (0°–90°N) with a latitudinal resolution of 5°. In each latitudinal belt, the surface is represented as a time-dependent ratio of seven different surface types (sectors): sea ice, ice-free ocean, ice-free continent, snowfield, North American ice cap, Eurasian ice cap, and Greenland's ice cap. The details of the model iterations are presented in section 2.1. This model was used successfully in the past in numerous studies investigating aspects of the glacial cycles, simulating ice volume over the last 500 kyr to 3 million years [Gallée et al., 1991, 1992; Berger and Loutre, 1997; Berger et al., 1998, 1999; Bar-Or et al., 2008]. At present, due to limited computing power, general circulation models (GCMs) cannot be used for continuous simulation of glacial cycles. Being simpler than a GCM, the LLN model includes the necessary components of the climate system needed to study glacial cycles, such as a description of ice sheet evolution over a million year time series, yet it is efficient enough to perform both sensitivity studies and forcing experiments over a large parameter space. 2.1. The Climate Model [7] Atmospheric dynamics in the climate model component of the LLN 2D are computed during each iteration using a dual layer quasi-geostrophic scheme. Precipitation and vertical and turbulent heat fluxes are taken into account as well. After a spin up of 15 years calculated with 1800 iterations of 3 days each, the climate model provides initial values (surface temperature, precipitation, etc.) for the ice sheet submodel. The ice sheet submodel in turn calculates the advance of the ice cap sectors using internal equations of motion, the equations of heat transfer, and the parameters provided by the climate model for initialization, in an online fashion updated at each iteration time step (1 kyr). In our experiments the total simulated time duration is set to 575 kyr, from 575 kyr before present (BP) until present time. 2.2. Modeling the Dust Radiative Effect [8] Since the radiation component of the LLN 2D assumes constant aerosol properties, we need to override the default aerosol values with time-dependent values corresponding to the historical dust flux record. For this purpose, we inspected the most detailed and continuous dust records available, taken from the ice cores drilled at the Vostok station and at the EPICA Dome C station in Antarctica [Jouzel et al., 1993; Petit et al., 1999; Delmonte et al., 2004]. Of the two records, we chose EPICA, because it extends as far back as 739 kyr BP. Although 232Th-derived marine sediment dust flux data from the equatorial Pacific are also available for the same periods, the marine sediment data lacks the required time resolution to simulate ice sheet processes [Winckler et al., 2008]. Furthermore, EPICA Community Members [2006] found a one-to-one coupling of glacial climate variability in Greenland and Antarctica, and Winckler et al. [2008] found a strong correlation between dust fluxes recorded in the tropics and in Antarctica. Winckler et al. [2008, p. 95] stated that, "the excellent correlation between dust fluxes recorded in the tropics and in Antarctica is particularly stunning given that the records from the two regions are based on different paleoarchives, and consequently independent age models, as well as different dust measures." Given the above and since there are no similar records available for the Northern Hemisphere, we utilize the same data (the EPICA Antarctica record) to force the both the Northern and Southern Hemisphere climates. (We do note that the amplitude of the dust record given by Winckler et al. [2008] is 2 orders of magnitude different from that of the EPICA record, and this is a source of uncertainty that cannot be resolved currently.) [9] In order to implement the variations in dust flux in our simulations, the atmospheric dust flux is assumed to be proportional to the dust accumulation data, using the age model and time sequence presented in the published data set [Delmonte et al., 2004]. We compute the simple ratio: Q(t) = F(t)/Fb, where F(t) is the dust accumulation record (in ppb according to the EDC2 time scale), and Fb is the background mean dust accumulation (91.15 ppb), represented by the mean EPICA dust accumulation from 1–17 kyr BP. Q(t) is therefore the nondimensional (relative) dust flux as a function of time. A similar unitless parameter of relative dust flux was used successfully in parameterizing time-dependent dust effects in the LLN model in our previous study [Bar-Or et al., 2008]. For our simulations including the direct radiative forcing of dust and the effect of dust on snow albedo in addition to the aerosol–Bénard cell effect, we set the optical thickness of continental type aerosols to be τ(t) = Q(t) · τc, where τc = 0.025 is the default (constant) optical depth of continental type aerosols in the model, and we set the albedo of snow and ice surfaces to be αs(t) = αs0 − 0.005 · Q(t), where αs0 = 0.85 is the default (constant) albedo of snow and ice surfaces in the model. 2.3. The Parameterization of Cloud Cover [10] As part of the radiation component of the LLN 2D climate model, global cloud cover is calculated every time step. Each latitudinal belt is represented by one effective cloud with averaged seasonal properties that remain constant during the run. Since there is no available paleodata of cloud cover, the latitudinal cloud fractions and properties are taken from past research, surveys, and simulations [Ohring and Adler, 1978; Berlyand and Strokina, 1980]. [11] Our modification of the model to incorporate the switching of opened Bénard cells above oceans into closed Bénard cells is achieved by the equation: where CF (j) is the new calculated cloud fraction for latitudinal belt j, CF0 (j) is the original model's cloud fraction for latitudinal belt j, B is the fraction of clouds above the ocean that are stratocumulus clouds, C is the relative cellular area that closes when the switch in direction of Bénard cell convection occurs, Q is the relative dust amount in the atmosphere compared to today's background, OF (j) is the oceanic surface fraction in latitudinal belt j as calculated within the model itself, and A(Q) (Figure 1) is the function that describes the frequency of the mechanism and is defined as follows: where Qc is the critical value of relative dust amount, which is the minimal value needed for the Bénard cell direction switch to occur, S is the exponential constant that controls the curvature of the exponential growth, and Θ(Q − Qc) is the Heaviside step function, defined as: Figure 1Open in figure viewerPowerPoint The intensity function A(Q) with inner parameters Qc = 1 and S = 1/2. [12] There are two main reasons for selecting this form of A(Q) over a simple step function that turns the mechanism on and off when the dust levels cross a certain critical value: (1) In our model, the ocean is homogeneous, and since the switching mechanism never works on 100% of the surface suitable for it to work on, we need to represent the ocean fraction that is affected by it; even in times of high dust concentrations in the atmosphere, there are some areas that stay clean and do not work according to the mechanism. (2) From recent observations [Rosenfeld et al., 2006; Wood and Hartmann, 2006; Wang and Feingold, 2009a, 2009b], the efficiency with which the Bénard cell switching mechanism occurs varies nonlinearly with aerosol concentration, increasing more rapidly for low aerosol concentrations than for high aerosol concentration. Such behavior is well described by the function A(Q) presented above (Figure 1), representing the relative dust concentrations. Typical values of the intensity function A(Q) can been seen in Figure 6. [13] Note that dust aerosols are not generally considered to be hygroscopic. However, recent studies have shown that while the activation of dust particles into cloud drops is slower than the activation of more soluble particles, such as ammonium sulfate, larger dust particles (greater than 2 μm in diameter) activate into cloud drops regardless of their composition, and smaller dust particles (with diameters between about 0.6 and 2 μm) activate into cloud drops if they contain even small amounts of K+, Mg2+, or Ca2+ compounds [Kelly et al., 2010; Kumar et al., 2010]. Therefore, dust aerosols may serve as cloud condensation nuclei. [14] An extensive set of sensitivity tests, described in section 4, led us to choose the following default values for B, C, Qc, and S: The chosen values of B and C show good agreement with previous studies. Wood and Hartmann [2006] estimate the cloud cover of low marine boundary layer clouds in low-cloud dominated regions of the eastern Pacific to be 50%, which serves as an upper bound on the value of B. They also present case studies with a decrease in cloud cover from 99% for closed cells to 32% for open cells. Similarly, Wang and Feingold [2009a, 2009b] model a cloud fraction of unity for closed cells and a cloud fraction of 0.20–0.35 for open cells, both of which are consistent with our chosen value of C. Therefore, the full equation used in this work to modify cloud cover is: [15] Note that while the nonlinear parameterization presented above is based on robust physical arguments, on past observations and simulations, and on careful sensitivity tests, as with any nonlinear mathematical description, it may produce some unwanted nonlinear phenomena. For example, the dust flux recorded at EPICA may already include a nonlinear amplification of the actual global dust source. Using this record as a linear proxy for the dust source function, which then in turn is nonlinearly amplified to parameterize the aerosol effect, may cause an unrealistic feedback and an overestimate of the aerosol–Bénard cell effect. [16] In addition, since the LLN 2D model only has two dynamical layers and only one effective/mean cloud layer, it is impossible to include the potential effect of variations in cloud height and thickness in the parameterization explicitly. Such effects and the possible semidirect effect of dust aerosols (in which heating by the aerosols, rather than their cloud drop nucleation properties, alters the cloud fraction) are included to a certain extent only implicitly by basing our parameterization on the statistics of stratocumulus cloud fields and their mean global coverage. The lack of an explicit description of changes in cloud height and thickness adds a measure of uncertainty to the parameterization. Nevertheless, this parameterization is the best available, being the only physical description of the aerosol–Bénard cell effect that is feasible for current climate models. 3. Results [17] All simulations presented here are compared to a reference simulation of the last 600 kyr that did not include any dust forcing. This exact simulation was also used as a reference by Loutre and Berger [2000]. 3.1. Single Forcing: The Aerosol–Bénard Cell Effect [18] The results of the simulation with the aerosol–Bénard cell effect in place are shown in Figure 2. It is clear that the aerosol–Bénard cell effect alone leads to a strong negative radiative forcing. With this forcing, the mean annual temperature during the glacial peaks decreases by up to 2.5°C in comparison to the reference simulation (solid line), and the maximal ice volume during the LGM reaches the value 51.86 × 106 km3. In addition, the ice volume values from this simulation are closer to those derived from variations in the ocean sea level [Rabineau et al., 2006, 2007; Charbit et al., 2007] and SPECMAP data set [Imbrie et al., 1989] than the values from the reference simulation (Figure 3, bottom). Figure 2Open in figure viewerPowerPoint Variations over the last six glacial-interglacial cycles of (top) the simulated Northern Hemisphere annual mean surface temperature and (middle) the Northern Hemisphere continental ice volume, with cloud fraction above oceans forced by dust (dashed line). The reference simulation (solid line) is not forced by any dust changes. (bottom) The δ18 O record, based on the SPECMAP data set [Imbrie et al., 1989], is presented for comparison. Figure 3Open in figure viewerPowerPoint Variations over the last six glacial-interglacial cycles of (top) the simulated Northern Hemisphere annual mean surface temperature and (middle) the Northern Hemisphere continental ice volume, with all three dust forcing mechanisms (dashed line). The reference simulation (solid line) is not forced by any dust changes. (bottom) The δ18 O record, based on the SPECMAP data set [Imbrie et al., 1989], is presented for comparison. 3.2. Combined Forcing: The Direct Effect, the Effect on Snow Albedo, and the Aerosol–Bénard Cell Effect [19] One of the reasons for the large uncertainty associated with the net effect of dust on the climate system is that different dust forcing mechanisms act in different directions [Forster et al., 2007]. In order to better understand the magnitude of the aerosol–Bénard cell forcing relative to the direct radiative effect of dust [Shell and Somerville, 2007; Alpert et al., 1998; Yu et al., 2006; Bar-Or et al., 2008] and to the effect of dust on snow albedo [Warren and Wiscombe, 1980; Peltier and Marshall, 1995; Calov et al., 2005; Motoyoshi et al., 2005; Krinner et al., 2006; Bar-Or et al., 2008], we run an additional simulation in which the three effects of dust, (namely, the direct effect of dust on the total aerosol optical depth, the effect of dust on the snow albedo, and the aerosol–Bénard cell effect) are taken into account simultaneously. The results are shown in Figure 3. [20] In this combined simulation, the Northern Hemisphere average surface temperature reaches very low values during glacial periods (4.8 ± 1.0°C lower during the LGM than the reference simulation) and reaches slightly higher values during interglacial periods (1.2 ± 1.2°C higher than the reference simulation). The global continental ice volume values, however, are almost identical to the reference run (differing by only 1–2% during the LGM). The low temperature values during the glacial periods may result from an overestimation of the global dust flux from Antarctic ice core data (refer to section 2.2). On the other hand, our results also point to the intriguing possibility that severe average surface temperature changes could have occurred while preserving the observed continental ice volume. 4. Sensitivity Tests [21] Sensitivity tests regarding key parameters governing the direct effect of dust and the dust-albedo effect were conducted by Bar-Or et al. [2008]. Here we present a brief discussion of the parameters that govern the aerosol–Bénard cell effect described by equations (1)–(3). The ranges of values we tested for the parameters B and C in equation (1) are 0.2 < B < 0.3 and 0.6 < C < 0.75, respectively (Figures 4 and 5), selected according to the estimations of previous research [Rosenfeld et al., 2006; D. Rosenfeld, personal communication, 2007]. The sensitivity tests presented in Figures 4 and 5 reveal that both the mean surface temperature and the NH ice volume are relatively insensitive to the values of B and C. Given this, we choose B = 0.3 and C = 0.7 as our default values for these parameters, which agree well with previous research [Wood and Hartmann, 2006; Wang and Feingold, 2009a, 2009b], as described in section 2.3. Figure 4Open in figure viewerPowerPoint Model sensitivity to parameter B: variations over the last six glacial-interglacial cycles of (top) the simulated Northern Hemisphere annual mean surface temperature and (bottom) the Northern Hemisphere continental ice volume with parameter B varying between B = 0.2 (solid line), B = 0.25 (dashed line), and B = 0.3 (dotted line) and with C = 0.6, S = 1/2, and Qc = 2/3. Figure 5Open in figure viewerPowerPoint Model sensitivity to parameter C: variations over the last six glacial-interglacial cycles of (top) the simulated Northern Hemisphere annual mean surface temperature and (bottom) the Northern Hemisphere continental ice volume with parameter C varying between C = 0.6 (solid line), C = 0.7 (dashed line), and C = 0.8 (dotted line) and with B = 0.3, S = 1/2, and Qc = 2/3. [22] Likewise, the A(Q) function (refer to Figure 1) contains two parameters, Qc and S, which control the curvature of the exponential growth. The critical dust concentration value Qc, which sets off the aerosol–Bénard cell effect, must be lower than unity, since we observe the aerosol–Bénard cell effect in today's background dust concentration (corresponding to Q = 1). We set Qc to be 2/3, a value lower than today's background but high enough to avoid setting off the aerosol–Bénard cell effect when the dust concentration is extremely low. We tested a wide range of values for the parameter S (0.1–2.0) (Figures 6 and 7). From Figure 6, there is a strong correlation between the intensity function and the dust flux record for all values of S tested. Figure 6Open in figure viewerPowerPoint Model sensitivity to parameter S: variations over the last six glacial-interglacial cycles of the normalized dust record (1/Qmax)Q (thick blue line) and the intensity function A(Q) for S = 0.1 (dotted black line), S = 0.5 (dashed black line), and S = 1 (solid black line). Figure 7Open in figure viewerPowerPoint Model sensitivity to parameter S: variations over the last six glacial-interglacial cycles of (top) the simulated Northern Hemisphere annual mean surface temperature and (bottom) the Northern Hemisphere continental ice volume with parameter S varying between S = 1/8 (dash-dotted line), S = 1/4 (dashed line), S = 1/2 (solid line), and S = 1 (dotted line) and with B = 0.3, C = 0.7, and Qc = 2/3. [23] In addition, from a separate set of sensitivity runs (Figure 7), we find that the model results, particularly the continental ice volume, are relatively insensitive to the value of S for values of S equal to or smaller than 1/2; that is, 1/2 is the highest value of S that will preserve the robustness of the model. Therefore, we choose S = 1/2 as the default value for this parameter in our model runs. [24] Note that the tests presented in Figures 4–7 show a higher sensitivity to changes in the tested parameters B, C, and S during glacial periods than during interglacial periods. This is consistent with the nature of the described mechanism, which operates only when the dust concentration exceeds a critical value (Qc) (which is more likely during glacial periods). 5. Discussion and Summary [25] In this work, we investigate the role of dust in glacial-interglacial cycles as a critical radiative forcing element using the LLN climate model. We modify the model to include an additional dust-dependent radiative forcing mechanism that has not been taken into account in this model so far, and has never been taken into account in multiglacial cycle simulations, namely the aerosol–Bénard cell effect. [26] With a time-dependent forcing, we find that the aerosol–Bénard cell effect plays a potentially important role in the glaciation mechanism, increasing the continental ice volume during the earlier glacial cycles to values closer to those simulated for the last glacial cycles. Combining the aerosol–Bénard cell effect with the direct radiative effect of dust and with the effect of dust on snow albedo reveals that when all three mechanisms operate simultaneously, the surface temperature variations between glacial and interglacial periods increase significantly and reach values of more than 10°C, almost twice the largest estimated difference in previous literature. However, the ice volume values remain almost identical to those in the reference run. This result points to the possibility that the natural surface temperature differences between glacial and interglacial periods were much larger than previously estimated, while ice volume values remain as observed from many natural sources (i.e., sea level differences, δ18 O records, etc.). [27] The purpose of the present, rather idealized, study is to highlight a single feedback rather than to attempt a completely realistic simulation. Being a nonlinear complex problem approached for the first time, this research still holds uncertainties. A higher temporal and vert