2.1 LOVECLIM

LOVECLIM (version 1.3) is a coupled Earth systems model of intermediate complexity (EMIC), which consists of atmosphere (ECBilt), ocean with dynamic sea ice (CLIO) and vegetation (VECODE) modules. It is fast enough to simulate the last glacial inception (120 to 100 ka) in less than 3 weeks using a single computer core. Therefore, it has been used to simulate a wide range of different climates from the last glacial maximum (Roche et al., 2007) through the Holocene (Renssen et al., 2009) and the last millennium (Goosse et al., 2005) to the future (Goosse et al., 2007).

Atmosphere

The atmospheric component (ECBilt, Opsteegh et al., 1998) is a spectral global quasi-geostrophic model, with T21 truncation, three vertical layers at 800, 500, and 200 hPa, and a time step of 4 h. The quasi-geostrophic structure of the model limits its ability to simulate equatorial variability and, hence, atmospheric interactions between the tropics and higher latitudes. To partially compensate for this, it has additional ageostrophic terms to improve the representation of Hadley cell dynamics (Opsteegh et al., 1998). Precipitation is computed from the precipitable water of the first layer according to a precipitation threshold for relative humidity (default 85%). The model contains simple schemes for short-wave and long-wave radiation, with radiative cloud cover prescribed by default (Haarsma et al., 1996).

Ocean

The oceanic component (CLIO – Coupled Large-scale Ice Ocean) is a 3-D primitive equation model with Boussinesq and hydrostatic approximations. The model is discretized horizontally on a 3^∘ × 3^∘ Arakawa B-grid, with 20 vertical levels on a z coordinate. This coarse-resolution enables CLIO to run fast enough for glacial cycle simulations. A free surface and a parameterization of down-sloping currents (Campin and Goosse, 1999) enables CLIO to receive freshwater fluxes and capture some of their impacts on dense water flows off continental shelves. Goosse et al. (2001) describe the model in detail. A major limitation of this model (and challenge for many GCMs) for paleoclimate studies is that the bathymetry and land mask can not be changed during a transient run (specifically, there is no available nor described implementation that can do so).

Sea ice

The sea ice component of CLIO is an updated version of the Fichefet and Maqueda (1997) dynamic–thermodynamic sea ice model. A visco-plastic rheology (Hibler, 1979) is used for horizontal stress balance. The thermodynamic component of the sea ice model considers sub-grid sea ice and snow cover thickness distribution, in addition to ice and snow sensible and latent heat storage.

Vegetation

VECODE is a dynamic terrestrial vegetation model with a simplified terrestrial carbon cycle (Brovkin et al., 2002). The model simulates the dynamics of two plant functional types (trees and grasses), in addition to deserts, and evolves their grid-cell fractions. These fractions are determined by the contemporaneous climate state and terrestrial carbon pool. More details about the model can be found in Brovkin et al. (1997, 2002).

LOVECLIM has been tested for both interglacial and glacial contexts. Nikolova et al. (2013) found that the large-scale changes in climate simulated by LOVECLIM for the last interglacial were in approximate agreement with those indicated by available proxies (differences in the ±2.5 ^∘C range for summer, and the −5 to 0 ^∘C range for winter). These changes were also similar to that from a full complexity atmosphere–ocean general circulation model (CCSM3). However, due to stronger polar amplification in LOVECLIM, smaller sea ice extent and higher surface temperatures are simulated in LOVECLIM compared to CCSM3 for the interglacial.

During the Last Glacial Maximum, LOVECLIM overestimates both the minimum and maximum Southern Hemisphere sea ice cover compared to paleo-proxy data, while CCSM only overestimates the minimum sea ice extent (Roche et al., 2012). Roche et al. (2007) also found a reasonable agreement between the atmospheric and oceanic estimates of LOVECLIM and proxy data during the LGM (e.g., disappearance of much of the PD Siberian boreal forest, seasonal sea-ice extent, sea surface temperature). However, the Atlantic deep ocean circulation is stronger in their simulation, which opposes the general inference of weaker AMOC during the LGM.

2.2 GSM

The GSM is built around a thermo-mechanically coupled ice sheet model. It includes a 4 km deep permafrost-resolving bed thermal model (Tarasov and Peltier, 2007), fast surface drainage and lake solver (Tarasov and Peltier, 2006), visco-elastic bedrock deformation (Tarasov and Peltier, 1997), positive degree day surface mass balance with temperature dependent degree-day coefficients derived from energy balance modelling results (Tarasov and Peltier, 2002), sub-grid ice flow and surface mass balance for grid cells with incomplete ice cover (Morzadec and Tarasov, 2015), and various ice calving schemes for both marine and proglacial lake contexts (Tarasov et al., 2012). For the results herein, ice shelves are treated using a crude shallow ice approximation with fast sliding. The GSM is run at a 0.5^∘ longitude by 0.25^∘ latitude grid resolution.

The GSM has three new features that have not been previously documented. First, ice calving has been upgraded to the more physically based scheme from DeConto and Pollard (2016). However, our implementation imposes the additional condition that the ice cliff failure mechanism is only imposed at ice marginal grid cells. Second, a temperature-dependent sub-shelf melt scheme that also depends on adjacent subglacial meltwater discharge from the grounded ice sheet has been added. The melt is proportional to the water temperature to the power 1.6 and to proximal subglacial meltwater discharge following the Greenland fjord modelling results of Xu et al. (2013). We also impose a quadratic dependence on ice thickness to concentrate sub-shelf melt near deep grounding lines in accord with the results of process modelling (e.g., Jacobs et al., 1992). Finally, a first-order approximation to geoidal deflection is now included. Details of these schemes will be in an upcoming submission that fully describes the revised GSM.

2.3 Model initialization

For the results herein, PD and glacial inception model runs are initiated with PD ice sheet thickness. The initial bed-thermal temperature field is set to the PD resultant field from a mix of best-fit past calibrated and ongoing calibration model runs (i.e., without LOVECLIM as in Tarasov et al., 2012 for North America).

To initialize the temperature field of existing ice sheets, results from previous transient runs are usually used in the GSM. However, this does not work for a cold start from PD fields, as the PD ice thickness fields will not line up fully with model results. As such, the ice temperature field is initially linearly interpolated from surface temperature to a basal temperature of −3 ^∘C. This enforces a frozen base to ensure a smooth spin-up and also provides warm enough ice to generate significant ice velocities. Ice velocity fields are then computed. The ice thermodynamics are subsequently partially spun-up over 5000 years (with fully coupled bed-thermal evolution). When the results of pre-Eemian Greenland and Antarctic calibrations become available, this will be used to initialize the respective Eemian ice sheet temperature and ice thickness fields.

Glacial inception surface topography is also offset from PD by the amount required to remove PD topographic discrepancies (to observed) from some past best fit calibration model runs.

As detailed below, the climate model is spun-up over an ensemble parameter dependent time interval prior to onset of the coupled model run.

3.1 Atmosphere to ice

At the end of each coupling time step, the coupler receives climate fields averaged over the last 10 years from ECBilt and converts them to monthly-mean values. These fields include the following:

• 2 m near surface air temperature and standard deviation;

• vertical 2 m temperature gradient;

• precipitation;

• evaporation;

• and latitudinal and longitudinal components of wind and the standard deviation of each.

LOVECLIM computes both 2 m near surface air temperatures (T2m) and surface (or skin) temperatures. We note at least one previous study has indicated the usage of LOVECLIM surface temperature for ice sheet modelling contexts (Roche et al., 2014), which we find problematic. Surface melt determination using positive degree days requires the T2m, and rain/snow fraction determination will be more accurately estimated with the T2m than the surface temperature. Ice thermodynamics would properly use the surface temperature, but this can be alleviated in part if the ice sheet model limits surface temperatures to 0 ^∘C over ice and snow for the ice thermodynamics.

Figure 2Total ice volume in sea level equivalent (m) at last glacial inception, synchronously coupled with 100-year, 20-year, and 10-year time steps.

Download

Given the simplified boundary layer physics of LOVECLIM, it may be that some weighted average of its T2m and surface temperature is a more appropriate estimate of “true” 2 m temperature. As shown in the Supplement, a raw average gives somewhat better overall fits to ERA40 2 m temperatures over Greenland and Antarctica but worse fits for July over North America and especially Eurasia using the default LOVECLIM tuning. Given these mixed results (and the possibility that after retuning the average of T2m and surface temperature would give better fits), we provide an option in the coupler to extract this average temperature from LOVECLIM instead of T2m.

Figure 3Vertical temperature lapse rate calculated by the coupler at PD over NA in (a) February and (b) July. (c) Shows the ice thickness difference between dynamic and constant 6.5 K km⁻¹ lapse rate (control) runs after running for 2 kyr, starting from the same 110 ka configuration. Black contours show the ice thickness in the control run. Thick black and green contours show the ice margin in the control and dynamic lapse rate run, respectively.

Download

The large difference in spatial resolution of the two models necessitates horizontal and vertical downscaling of the climatic fields. The GSM receives climatic fields on the LOVECLIM grid, and downscales them to its own grid resolution using bilinear interpolation.

The downscaled standard deviation of temperature (using 4-hourly ECBilt data for each month averaged over the last 10 years of each coupling time step) is used to compute monthly positive degree days, with the usual assumption of a Gaussian distribution around the monthly mean. This is opposed to the traditional practice of assuming a constant value, usually between 5 and 7 ^∘C.

Figure 4Impact of advective precipitation downscaling inclusion in the coupled model; NA ice thickness difference at 110 ka between simulations with and without the advective precipitation method. Contours show the ice thickness in the control run. Thick black and green contours show the ice margin in the control and advective precipitation run, respectively. LOVECLIM parameters are set to default values.

Download

Figure 5The upstream ocean temperature profile sites and corresponding downstream sectors assigned to these profiles for ocean–ice coupling in (a) NA and Greenland and (b) EA.

Download

3.1.1 Vertical temperature gradient

Large grid resolution differences between ECBilt and the GSM result in surface elevation differences between the two models, especially in places with steep topography. The altitude dependence of temperature in such regions can drastically affect the type of precipitation and surface mass balance of the ice sheet. Therefore, in addition to horizontally downscaling the temperature from LOVECLIM to the GSM, a vertical correction of temperature is required.

By monitoring 25 sites spread over a 15 650 km² area and with an altitude range from 130 to 2010 m on the Prince of Wales Icefield for 2 years, Marshall et al. (2007) found a mean daily vertical surface temperature gradient of −4.1 K km⁻¹, with an average summer gradient of −4.3 K km⁻¹. These values are less than the standard mean free-air temperature lapse rate that is often used for extrapolations of sea level temperature to higher altitudes (−6.5 K km⁻¹) (e.g., Glover, 1999; Arnold et al., 2006; Raper and Braithwaite, 2006). Marshall et al. (2007) also find a vertical surface temperature gradient of −6 to −7 K km⁻¹ on steep regions in summer, and around −2 K km⁻¹ in regions where northerly anticyclonic flow is more common. In addition, Gardner et al. (2009) find significant spatio-temporal variations in vertical temperature gradients across four glaciers in the Canadian High Arctic.

The GSM uses the near-surface vertical T2m gradient calculated by the coupler at the end of each time step to downscale the temperature field over its high-resolution grid. In each LOVECLIM grid cell, the coupler first determines the highest and lowest elevations from the GSM topography constrained by the cell's boundary. Next, the T2m for these two elevations is calculated using the inherited scheme from the LOVECLIM atmospheric model (as detailed in Roche et al., 2014). The resulting temperature and elevation difference between the two points is then used to calculate the temperature lapse rate in that LOVECLIM grid cell.

Figure 6Ice thickness difference from the control run (dynamic ocean temperature) at 108 ka for three test runs starting from the same 110 ka control run restart. The three test runs are as follows: (a) PD ocean temperature run, (b) the fixed ocean temperature at −2 ^∘C run, and (c) the temperature averaged over ocean layers run. Contours show the ice thickness in the control run. Thick black and green contours show the ice margin in the control and the other run in each of the panels, respectively.

Download

Figure 3a and b show the present-day vertical T2m lapse rate calculated by the coupler for summer and winter. The derived lapse rate has strong spatial and temporal variation over NA and Greenland. The impact of this variation is shown in Fig. 3c. Starting from the same 110 ka configuration, the difference in ice thickness after 2 kyr between a dynamic temperature lapse rate run and a control run (default LOVECLIM parameters) with a 6.5 K km⁻¹ lapse rate can reach over 1 km.

Evaluating the appropriateness of our vertical temperature downscaling approach is difficult, especially when considering glacial/interglacial changes. Using a global climate model (CCSM3), Erokhina et al. (2017) found significantly larger surface slope lapse rate values over the Greenland ice sheet during the LGM compared to pre-industrial values (February mean increase of about 3.7 K km⁻¹ and about 0.9 K km⁻¹ for July). In contrast, our T2m mean LGM lapse rate over Greenland is 0.8 K km⁻¹ stronger for February and 0.2 K km⁻¹ weaker for July compared to that of the PD. However, neither lapse rate is an accurate a priori choice for vertical downscaling. A need remains for a multi-resolution modelling study to compare a “true” downscaling vertical temperature gradient with the various possible lapse rates that can be derived from a single resolution model.

Figure 7The impact of meltwater runoff inclusion in the coupled model; (a) total ice volume evolution at glacial inception with (green) and without (red) dynamic meltwater routing, and (b) NA ice thickness difference at 110 ka with and without dynamic runoff routing. Contours show the ice thickness in the simulation without dynamic runoff routing. Thick black and green contours show the ice margin in the control and dynamic runoff routing run, respectively.

Download

3.1.3 No bias correction

Studies from Goosse et al. (2007), Mairesse et al. (2013), Renssen et al. (2009), Widmann et al. (2010), Roche et al. (2007), van Meerbeeck et al. (2009), and Otto-Bliesner and Brady (2010) demonstrate LOVECLIM's overall ability to simulate the last millennium, the Holocene, and the Last Glacial Maximum (LGM) climates in agreement with observed and proxy records. However, the model still suffers from a high temperature bias at low latitudes, an overly symmetric distribution of precipitation between the two hemispheres, an overestimation of precipitation and vegetation cover in the subtropics, weak atmospheric circulation, and an overestimation of the ocean heat uptake over the last decades (Goosse et al., 2010).

Table 2Ensemble parameters that are varied in the historical transient simulation ensemble. Column 2: the distribution of parameter values versus their range for each parameter. Column 3: change in 1950–1980 CE mean summer 2 m temperature and winter precipitation over four selected regions when each parameter is varied independently from its minimum to its maximum value. In each sensitivity run, all the other parameters are fixed to LOVECLIM default values, with a spin-up length of 4000 years, simple upscaling method, default LOVECLIM cloud radiative forcing, a start year of 1500 CE, and a dynamic vertical temperature lapse rate.

Download Print Version

The extent to which these biases are due to the tuning of LOVECLIM parameters and missing couplings with the rest of the Earth/climate system is unclear. Therefore, we do not apply a bias correction to atmospheric fields and instead examine the extent to which an ensemble parameter sweep can reduce the bias. As detailed below, a reduction in the PD regional temperature and precipitation bias occurs for various members of our perturbed parameter ensemble.

The control run (with all LOVECLIM parameters set to their default values, and other coupling parameters as described in the caption of Table 2) shows the highest temperature bias in the “South NA” region (∼5 ^∘C), with slightly colder temperatures in the “North NA” (∼1 ^∘C). The temperature bias over EA is less significant, and is also less latitude dependent (both “North EA” and “South EA” are biased by <2 ^∘C). A reduction in the regional temperature and precipitation bias is observed in various members of our ensemble of simulations for the PD, as introduced below. The regional temperature and precipitation bias relative to observed (Table 3) over NA and EA can reach zero for some ensemble members for both summer and winter. Although there is no individual run with zero bias in all the regions, a number of selected runs show reduced temperature biases (between −1 and 1 ^∘C) in all four regions compared to that of the control run.

Table 3The sieved sub-ensemble and observed mean summer and winter 2 m temperature and precipitation averaged over four latitudinal bands for the 1950–1980 CE interval.

Download Print Version | Download XLSX

3.2 Ocean to ice: sub-shelf melt

Sub-ice shelf melt is a challenge for paleo coupled ice sheet climate modelling given the dependence on unresolved basin-scale circulation. As a first-order approximation, we assume that upstream ocean temperature at the same depth corresponds to the local sub-shelf temperature. To facilitate fast and simplified coupling, given the complexity of ocean grids in most ocean general circulation models, we only extract upstream ocean temperature vertical profiles from LOVECLIM at the end of each coupling time step for a number of chosen index sites as indicated in Fig. 5 and use these for downstream marine sectors. We selected these sites (seven over NA + Greenland and four over EA) by examining the PD ocean temperature climatologies from CLIO (at various depths) while taking ocean currents into account. Our site selection was predicated on the constant bathymetry and land mask of CLIO and would need updating for a model with dynamic land mask/bathymetry. The downstream masks for the profile sites extend onto land where applicable when grounding line retreat beyond the fixed ocean mask of CLIO (i.e., onto the land mask) is possible.

To test the impact of this regional disaggregation of ocean temperatures, we generated three test cases: ocean temperature forcing set to PD value, to −2 ^∘C, and to the contemporaneous average across the above index sites. Starting from a 110 ka restart, all three options have local ice thickness differences greater than 1 km after 2 kyr compared to that with the standard coupling (Fig. 6).

3.3 Ice to atmosphere

Changes in both the topography and the ice mask can affect the global circulation patterns by influencing the stationary waves and the jet stream. At the end of each coupling time step, the coupler receives the updated topography and ice thickness fields from the GSM. The topography field is upscaled to the ECBilt grid and then used for the next LOVECLIM run step. Given the large difference in grid resolution, the choice of upscaling scheme is not clear a priori. Therefore, we have implemented three different schemes to upscale the topography from the GSM high-resolution grid to the ECBilt low-resolution grid.

3.4 Ice to ocean

3.4.1 Topographically self-consistent and mass conserving freshwater discharge

The melting of continental ice sheets provides a freshwater source to the ocean that affects global sea level and the AMOC. Dynamical ocean models indicate that the strength of the AMOC in the North Atlantic is sensitive to the freshwater budget at the sites of the formation of North Atlantic Deep Water (Rahmstorf, 1995).

Figure 8Selected north and south zones over North America and Eurasia for PD sensitivity analysis.

Download

As the GSM self-consistently computes surface drainage (while conserving mass) for the evolving topography while LOVECLIM surface drainage is hard-coded for PD topography, precipitation within LOVECLIM is masked out where covered by the GSM grid. The coupler then distributes GSM freshwater ocean discharge to corresponding LOVECLIM ocean discharge grid cells. In areas where the GSM grid edge is terrestrial, the GSM discharge is added to the corresponding LOVECLIM grid cell for internal runoff routing. As topographic gradient changes from glacial isostatic adjustment (GIA) are small outside of the GSM grid, this scheme should give results close to what would be achieved with a drainage solver using a topography globally subject to GIA.

We describe the runoff routing in the coupler and the connection between the GSM and LOVECLIM drainage basins in more detail in the Supplement.

Figure 7a provides an example of how the inclusion of meltwater runoff in the coupled model improves ice sheet growth at glacial inception. Although the impact is small at the first stages of ice formation due to small ice volumes with negligible runoff rate changes, ice in the simulation including runoff grows faster as it gains more volume. The difference reaches its maximum at 110 ka, with about 50 % more ice in the run with dynamic drainage routing, including much thicker and more extensive ice over North America (Fig. 7b). Compared to the control run (no dynamic drainage routing), the AMOC strength drops by 15 %, and the sea ice extent shows an increase of 5 % in winter and ∼15 % in summer (not shown). The combination of these AMOC and sea ice changes yields a cooler summer in the run with dynamic drainage routing, which results in less ice sheet melt.

3.4.2 Bering Strait

The Bering Strait is a narrow strait with a present depth of approximately 50 m between Siberia and Alaska, through which relatively fresh North Pacific water is transported to the Arctic. From there, the North Pacific water is transported to the Greenland Sea and the North Atlantic. This less-saline water affects the upper ocean stratification and thus the strength of deep ocean convection and the AMOC, which in turn, has impacts on the global climate (Shaffer and Bendtsen, 1994; De Boer and Nof, 2004; Hu et al., 2008).

Due to the ice sheet growth and associated sea level lowering, the Bering Strait was often closed during glacial cycles, limiting the freshwater flow from the Pacific Ocean to the Arctic. LOVECLIM does not explicitly compute the direct connection between the Pacific and the Arctic through the Bering Strait, so the transport is parameterized by a linear function of the cross-strait sea level difference in accordance with geostrophic control theory (Goosse et al., 1997). The coupler interpolates the Bering Strait scaling at each coupling step between the PD value (0.3, 50 m depth) and a closed strait (0.0, 0 m depth) using the relative sea level at the Bering Strait as computed by the GSM. Given the shallowness of the strait, the accuracy of the GSM in representing sea level changes (based on its viscoelastic bedrock response and first-order Geoidal correction) plays a potentially important role here.

Figure 9Total ice volume at the last glacial inception with the cloud parameterization (green) and the PD cloud cover forcing (red).

Download

4.1 Snow and ice albedo

Changes in the snow and ice area and type have an amplifying effect on climate by modifying the surface albedo. During summer, the balance between absorbed and reflected solar energy at the ice sheet surface is the dominant factor controlling surface melt variability in the ablation zone (van den Broeke et al., 2008). The parameterization of the surface albedo in LOVECLIM takes the state of the surface (frozen or melting) and the thickness of the snow and ice covers into account (Goosse et al., 2010). We include all types of snow and ice albedo (i.e., snow, melting snow, and bare ice) in our ensemble parameter set.

The snow albedo, bare ice albedo, and melting ice albedo rows in Table 2 show the range of albedo values for each type and their climate sensitivity. Increasing snow albedo results in a reduction of winter precipitation over all four regions with an extended effect over summer temperatures, as expected. However, the same feedback is not as straightforward for bare ice albedo and melting ice albedo. Although the increase in bare ice albedo shows an expected cooling effect over all regions with a smaller influence on precipitation, an increase in melting ice albedo causes all regions to get between 4 and 6 ^∘C warmer and increases the winter precipitation.

4.2 Climate initialization/spin-up

Before starting a coupled transient climate simulation, it is necessary to allow the atmosphere and ocean to adjust to the initial boundary conditions and external forcings. Model spin-up, and therefore the initial state of the climate system, can be a major source of uncertainty in climate modelling especially given the millennial timescale of deep ocean circulation.

The general approach to spin-up the ocean is to run the ocean to an equilibrium state under fixed external forcings (e.g., Manabe et al., 1991; Johns et al., 1997). However, as the climate system is unlikely to ever be in equilibrium, this choice lacks justification. We include two parameters to control the initial state of the system: LOVECLIM spin-up start year, and LOVECLIM spin-up length. All spin-ups are performed using transient orbital and CO₂ forcings ranging from 3000 to 5000 years but without the GSM coupling. The combination of these two spin-up control parameters results in slightly different coupled transient start times, each with different initial ocean and atmosphere states. For the runs herein, we constrain the spin-up to end between 1400 and 1600 CE. Increasing the spin-up length has a cooling and drying effect in the coupled model with PD boundary and initial conditions, while starting the transient coupled run from earlier years results in slightly warmer and wetter conditions (Table 2).

Figure 11Greenland ice thickness ensemble (a) mean, and (b) standard deviation at PD.

Download

4.3 Upscaling

The three different upscaling methods described in Sect. 3.3.1 are evenly distributed between ensemble members as shown in Table 2. By switching between three methods, we calculate the highest temperature and precipitation changes over four regions and plot them in the last column of Table 2. The highest temperature sensitivity to the upscaling method is recorded in NorthNA followed by the NorthEA zones, and the highest precipitation sensitivity is seen in the SouthEA zone.

4.4 Precipitation threshold

ECBilt only accounts for humidity, and thus precipitable water, between the surface and the 500 hPa layer. Above 500 hPa, the atmosphere is assumed to be dry, meaning that all the water transported by atmospheric flows into this region precipitates. Below the 500 hPa layer, ECBilt precipitates all the excess water above a fixed threshold (default 0.83) multiplied by the vertically integrated saturation specific humidity (Goosse et al., 2010). This parameter has the largest relative impact for NorthEA temperature (4 ^∘C over the parameter range equivalent to 60 % of mean).

4.5 Cloud radiation parameterization

The representation of clouds is one of the largest sources of uncertainty in models. They play an important role in regulating the surface energy balance of ice sheets, with competing warming and cooling effects at the surface through changes to short-wave and long-wave radiative fluxes. The effect of ice sheets on cloud formation is also significant. The growth of ice sheets results in tropospheric cooling and a reduction in humidity. This colder and drier troposphere displaces the upper tropospheric stratiform clouds downward, and reduces the low level stratiform cloud cover around the ice sheets (Hewitt and Mitchell, 1997).

The total downward and upward long-wave radiative scheme in ECBilt is a function of the vertical profile of the temperature, the concentration of various GHGs, and the humidity, and is computed for both clear sky and cloudy conditions. The radiation computed for each grid cell is then the weighted average of these two conditions based on the cloud coverage. The default ECBilt configuration prescribes radiative cloud coverage to the PD ISCCP D2 dataset (Rossow, 1996). The total downward and upward short-wave radiative fluxes depend on the transmissivity of the atmosphere, which also relies on the prescribed cloud cover (Goosse et al., 2010).

Figure 12The distribution of global annual mean 2 m temperature difference between ensemble members and observations averaged from 1950 to 1980 CE. Grey and green bars represent the 2000 ensemble members and the top-performing 500 ensemble members, respectively.

Download

Given the importance of cloud radiative feedbacks on ice sheet evolution, the use of a prescribed PD cloud cover for paleoclimate modelling lacks justification. Therefore, we have added a simple cloud parameterization scheme similar to the precipitation parameterization scheme described in Sect. 4.4. The only difference here is the humidity threshold for cloud formation, which is assumed to be 10 % less than the precipitation threshold, allowing cloud cover without precipitation. Including the dynamic cloud cover radiation feedback in the coupled model slightly decreases the total ice volume at glacial inception (Fig. 9) through reduced humidity during glacial conditions reducing the cloud cover.

As evident in the “Temperature and precipitation sensitivity” column in Table 2, regional temperature and precipitation are sensitive to each of our ensemble parameters in the coupled model. However, due to the non-linearity of the climate system, the combined effect can be significantly different. In the next section, we explore the coupled model response to all these parameters in an ensemble of simulations.

5.1 The 2 m temperature and precipitation

The ensemble distribution of the annual mean global T2m anomaly with respect to observations shows that the majority of the ensemble members fall within ±3 ^∘C of observations (grey bars in Fig. 12). However, as ice sheet build up is a function of both temperature and precipitation, most of the warm biased simulations fail to maintain ice-free conditions over NA and EA during PD due to a high winter precipitation bias (Table 3).

Our four latitudinal bands defined in Sect. 4 (Table 3) provide more relevant temperature metrics for NH ice sheet contexts. All four regions have higher ensemble mean seasonal T2m and precipitation compared to observations. However, the observations are covered well within two standard deviations for all regions, and temperature is covered within one standard deviation for most regions.

The seasonal cycle provides a partial test of a model's response to orbital forcing on Milankovitch scales. The ensemble mean seasonal cycle (difference between mean summer and mean winter) is within one standard deviation of the reanalysis data for all regions and for both temperature and precipitation (Table 3). Furthermore, aside from NorthEA, the diagnostic regions have a mean difference between summer and winter ensemble temperatures within a degree Celsius of that of the reanalysis climatology.

5.2 Northern Hemisphere jet stream

Jet stream latitude and oscillations have a strong control over storm tracks and the boundary between polar and subtropical air masses. Therefore, they are critical factors in controlling where and when an ice sheet margin advances or retreats. Due to the low vertical resolution of ECBilt, we compare the ensemble zonal mean of the 200 hPa zonal wind (as opposed to the more usual 300 hPa diagnostic level) with observations in winter and summer over NA and EA (Fig. 13a). The ensemble shows good agreement in capturing the maximum zonal velocity, but there is a 10 to 15^∘ shift northward in the latitude of the jet in both seasons. This is likely due to the reduced temperature gradient between low and high latitude.

We also compare the 30–80^∘ N meridionally averaged meridional wind at 200 hPa of the ensemble mean and the observations to diagnose the Rossby waves. In the summer, the longitudes of the troughs and ridges from the Pacific Ocean to the Atlantic Ocean largely match the reanalysis output within ensemble range (red line in Fig. 13b) with the largest discrepancies over the Eurasian region. During the winter, although the general pattern of the jet stream oscillations still agrees between the model and the observations (troughs over Eurasia and North America), the mismatch between ensemble members and observations becomes more significant given the higher Rossby wave number of the model ensemble.

5.3 Sea ice

High latitude sea ice acts as an insulator for both heat and moisture between the atmosphere and ocean, which are the two controlling factors for terrestrial ice sheet surface mass balance. We use the area and minimum latitude extent of the NH sea ice as relevant diagnostics. The general warm bias of the ensemble is reflected in the reduced ice area of the ensemble for both seasons, barely capturing the observed area within the one standard deviation range of the ensemble (Fig. 14). Both March (maximum) and September (minimum) NH sea ice areas show gradual decreases in the ensemble mean as the greenhouse gas concentration increases in the model.

Our filter condition for our sub-ensemble still permits a wide response of modelled components. For instance, averaging from 1950 to 1980 CE, the Pacific Ocean sea ice shows higher sensitivity to ensemble parameters during its maximum seasonal extent than the Atlantic. The sea ice minimum latitude in the Pacific ranges from 60 to 45^∘ N (not shown), in comparison to the observed value of 60^∘ N (Walsh et al., 2015).

5.4 AMOC

The AMOC transports large amounts of heat and salinity between high and low latitudes. Both paleoclimate proxy records (McManus et al., 2004) and climate model simulations (Liu et al., 2009) show that the AMOC experiences significant changes over a glacial cycle.

Our “bounding reality” criteria is not met for at least two AMOC features of the sub-ensemble. The PD ensemble mean AMOC strength is weaker than the reanalysis data from European Centre of Medium-Range Weather Forecasts (ECMWF ORA-S3, Balmaseda et al., 2008). The temporal mean of the reanalysis data is only captured by the maximum ensemble range (Fig. 15a). The ensemble mean shows a slight increase in the AMOC strength from 1965 to 1980 CE, which is not seen in ORA-S3. Furthermore, the temporal variability of the CLIO AMOC lacks the strong amplitude of the low frequency component of observations (as displayed by the maximum and minimum – time averaged – AMOC strength runs in Fig. 15a).

The maximum AMOC stream function strength is seen around 50^∘ N at 1 km depth. The ensemble variance is also highest in the same region, in addition to 0^∘ latitude at the same depth (Fig. 15b).