Coupling the Glacial Systems Model ( GSM ) to LOVECLIM : description , sensitivities , and validation

We have coupled an Earth Systems Model of Intermediate Complexity (LOVECLIM) to the Glacial Systems Model (GSM). This coupling includes a number of interactions between ice sheets and climate that are often ignored: dynamic meltwater runoff routing, novel down-scaling for precipitation that corrects orographic forcing to the higher resolution ice sheet grid (“advective precipitation”), dynamic vertical temperature gradient, and ocean temperatures for sub-shelf melt. The 5 sensitivity of the coupled model with respect to the selected parameterizations and coupling schemes is investigated. Each new coupling feature has a significant impact on ice sheet evolution. An ensemble of runs is used to explore the behaviour of the coupled model over a set of 2000 parameter vectors using Present-Day (PD) initial and boundary conditions. The ensemble of coupled model runs is compared against PD reanalysis data for atmosphere (surface temperature, precipitation, jet-stream and Rossby number of jet), ocean (sea ice, sea surface 10 temperature, and AMOC), and Northern Hemisphere ice sheet thickness and extent. The parameter vectors are then narrowed by rejecting model runs (1700 CE to present) with regional land ice volume changes beyond an acceptance range. The selected sub-set forms the basis for ongoing work to explore the spatial-temporal phase space of the last two glacial cycles.


Introduction
Transitions between glacial and interglacial states have been a periodic feature of the Earth's climate for the last few million years.The driver of these transitions is understood to be orbital forcing (Berger, 2014;Birch et al., 2016;Birch et al., 2017;Hahn et al., 2015;Rind et al., 1989), with an important role for CO 2 variations (Elison Timm et al., 2015;Ganopolski et al., 2016)).Nevertheless, the role that climate feedbacks play in amplifying or inhibiting the responses to these forcings is not clear.
Given available proxy data, how well do we know the progression of these glacial cycles?Is there more than one way each transition could have occurred?How sensitive were these glacial cycles to small perturbations in the external forcings (e.g., volcanic eruptions)?In order to address these questions, we need to understand the relative importance of different feedbacks between ice sheets and other aspects of the climate system.We can build such understanding by probing this phase space with physically-based models that include the pertinent feedbacks on glacial timescales.
Temperature and net precipitation (the solid/liquid fraction thereof) encompass the main atmospheric impacts on ice sheets.
Marginal ice sheet surface mass-balance is very sensitive to the vertical temperature gradient.As indicated by the observations 2.An advective precipitation down-scaling scheme which accounts for wind direction and topographic slopes.
Table 1 compares the feedbacks included in previous coupled modelling studies to this one.There are two main feedbacks yet to be implemented.First, the dust cycle and its impact on atmospheric radiative balance and ice surface albedo (and therefore surface mass balance) awaits future work.Second, the LOVECLIM ocean component can not handle changing bathymetry and landmask.It does have a parameterized Bering Strait throughflow which permits shutdown of throughflow when local water depth approaches zero.
The climate model used in the coupled model needs to be fast enough to simulate tens of thousands of years in a reasonable time interval, while sufficiently complex to include all important climate dynamics.We tested every available fast model that included ocean, atmosphere and dynamical sea ice components, and found a number of published models to be numerically unstable or otherwise unable to run or port.The only stable model with all these components was LOVECLIM.The paper is structured as follows.We first introduce the models in section 2. Next, we describe the coupling schemes between the ice sheet model and the atmosphere and the ocean models in section 3 and their impact on ice sheet evolution during the last glacial inception.In section 4, we introduce the ensemble parameters for the coupled model and the ranges used for each.We examine the sensitivity of the coupled model to changes in each parameter for Present-Day (PD) climate.Then we test the ensemble parameter set using our coupled model with PD initial and boundary conditions and compare the results with observational/reanalyzed data.
Ocean The oceanic component (CLIO: Coupled Large-scale Ice Ocean) is a 3D 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) describes the model in detail.The main downside of this model for paleoclimate studies is its inability to change bathymetry and land mask.

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, and ice and snow sensible and latent heat storage.
Vegetation VECODE is a dynamic terrestrial vegetation model with 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) and Brovkin et al. (2002).Nikolova et al. (2013) investigated the performance of LOVECLIM and a full-complexity atmosphere-ocean general circulation model (CCSM3) during the last interglacial, and found both models qualitatively simulate the large-scale climate changes in line with proxy records.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.

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, 2017), and various ice calving schemes for both marine and pro-glacial 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 0.5 o longitude by 0.25 o latitude grid resolution.
The GSM has three new features that haven't been previously documented.First, ice calving has been upgraded to the more physically based scheme of 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 sub-glacial meltwater discharge from the grounded ice sheet has been added.The melt is proportional to the water temperature and proximal sub-glacial meltwater discharge.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 fully describing the revised GSM.

Coupling
The coupler is designed to regrid and exchange data between the ice sheet model (GSM) and LOVECLIM (ECBilt and CLIO) in both directions with minimal adjustment to the model code.Figure 1 displays all the fields the coupler transfers between different component models and the processes involved.
Due to the computational costs of coupling and trivial variations of ice sheets over small time-scales, the ice model and the climate model run for a certain number of years before receiving updated fields from the other model.On the other hand, using large coupling time-steps can also introduce errors into the results.To test the effect of the coupling time-step on ice sheet evolution, we used three different coupling steps (100, 20, and 10 years) to simulate the last glacial inception starting at 120 ka.With similar boundary and initial conditions for all three simulations, the total ice volume growth rate shows very small variance (max 5%) between all runs before the first total ice volume peak at around 110 ka (Fig. 2).The 100-year coupling-step run (red line in Fig. 2), however, diverges from the other two during the retreat phase.Therefore, we choose 20-years as the optimum coupling step for all of our ensemble simulations.
The ice sheet model exchanges data with both the atmosphere and the ocean models at the end of each coupling step.The fields that are passed are described in detail in the following.-Latitudinal and longitudinal components of wind standard deviation 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 bi-linear interpolation.
The downscaled mean standard deviation of temperature (within a month) is used to compute Positive Degree Days for each month, with the usual assumption of a Gaussian distribution around the month mean.This is opposed to the traditional practice of assuming a constant value of about 5 • C.

Vertical temperature gradient
Large resolution differences between ECBilt and the GSM result in height 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 GSM, a vertical correction of temperature is needed.
By monitoring 25 sites spread over a 15,650 km 2 area and with an altitude range of 130 to 2010 m on the Prince of Wales Icefield for two years, Marshall et al. (2007) found a mean daily vertical temperature gradient of -4.1 • C/km, with an average summer gradient of -4.3 • C/km.These values are less than the moist adiabatic free-air temperature lapse-rate that is mostly used for extrapolations of sea level temperature to higher altitudes (-6.5 • C/km) (Glover, 1999;Flowers and Clarke, 2002;Geosci. Model Dev. Discuss., Arnold et al., 2006;Bassford et al., 2006b, a;De Woul et al., 2006;Raper and Braithwaite, 2006).Marshall et al. (2007) also find a vertical temperature gradient of -6 to -7 • C/km on steep regions in summer, and around -2 • C/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 temperature gradient calculated by the coupler at the end of each time-step to downscale the temperature field over its high resolution grid.To extract this near-surface vertical temperature gradient, we apply the same approach ECBilt uses to calculate the surface temperature, but for two virtual surfaces corresponding respectively to the lowest and highest GSM altitudes within the same grid cell, as suggested by Roche et al. (2014).ECBilt computes vertical temperature profiles based on reference temperature profiles for 19 virtual levels, 27 different regions, and 12 months from NCEP-NCAR reanalysis data (Kalnay et al., 1996), assuming temperature anomalies with respect to these temperature profiles vary linearly with the log of the pressure.

Advective precipitation downscaling
LOVECLIM calculates rain and snow for each grid cell based on its temperature fields.Since the GSM is using a different vertical temperature gradient as described in the previous section, the calculation of snow versus rain precipitation from LOVE-CLIM is not consistent with the climate on the GSM grid.A common approach is to linearly interpolate both precipitation and evaporation fields onto the high resolution ice sheet model grid, calculate the net precipitation amount, and finally determine the amount of rain and snow for each grid cell using the downscaled temperature.Here, we apply a new approach to precipitation downscaling to also account for orographic forcing at the ice sheet grid resolution.
The scheme assumes that orographic precipitation effects for upslope winds will be proportional to the vertical velocity induced by the surface slope and therefore to the dot product of the horizontal wind speed and surface slope (S GSM and S AT M ) with the latter given by: The k in the above equations indexes a representive range of wind vectors (u(x, y, month, k)) for each month.To simplify the coupling and still capture wind variation, we use monthly climatologies of mean wind velocity and its standard deviation in the determination of S GSM and S AT M .We compute the advective precipitation correction factor (f k p ) using in turn the S's as a function of: mean, and mean ± one standard deviation and then sum over these factors with appropriate weights (W (k)) for a Gaussian distribution.This correction is based either on the ratio of the S terms for S AT M > 0 (i.e., upslope winds) or else their difference (to transition into precipitation-shadowing).In detail, with the inclusion of a regularization term (µ, that governs the transition to precipitation-shadowing) and bounds (f pmin and f pmax ), this takes the form: The loop is carried out for each point on the GSM grid.The net correction for each corresponding point on the lower resolution atmospheric grid is then accumulated to generate a rescaling coefficient that is mapped back to each f p (x, y, month) on the GSM grid to ensure mass conservation.The scheme is currently implemented with µ = 0.005 and f pmin = 0.2 and f pmax = 5.0.
The new advective precipitation downscaling results in increased ice sheet volume and southern extent for the North American ice sheet during the inception phase.This increase is largest for the southeastern sector of the ice sheet (Fig. 3).Ice thickness also decreases in some regions due to precipitation-shadowing.  the two hemispheres, an overestimation of precipitation and vegetation cover in the subtropics, a weak atmospheric circulation, and an overestimation of the ocean heat uptake over the last decades (Goosse et al., 2010).
It is unclear 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.We therefore do not apply a bias correction to atmospheric fields and instead examine the extent to which an ensemble parameter sweep can reduce the bias.

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 basinscale 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 for a number of chosen index sites as indicated in Fig. 4 and use these for downstream ocean bands.
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 o C, and set to the contemporaneous average across the above index sites.Starting from a 110 kyr restart, all three options had local ice thickness differences greater than 1 km after 2 kyr compared to that with the standard coupling (not shown).

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 into the ECBilt grid resolution and then used for the next LOVECLIM run step.Given the large difference in grid resolution, the choice of upscaling scheme is not a priori clear.We have therefore implemented three different schemes to upscale the topography from the GSM high resolution grid into the ECBilt low resolution grid.

Topography Upscaling and Ice-mask
Simple average method In this method, the coupler simply calculates a weight for each high-resolution grid cell based on the fraction of the cell located inside the coarse grid cell.These weights are then used to calculate the average altitude of each ECBilt cell from GSM orography.
Envelope Method In the envelope method, a weighted standard deviation of the altitude of all GSM cells inside the ECBilt cell is added to the simple average altitude from the previous method.The envelope method works reasonably well to preserve the overall topographic peaks, but it can introduce a phase shift in the terrain field, broaden ridges, and raise the height of even relatively broad valleys.Here, H i,j is the model terrain height, ω is a predefined weighting factor (in our experiments 0.5), and σ i,j is the standard deviation at the model grid point.

Silhouette method
The silhouette method combines the simple average altitude with a silhouette height.The silhouette height is defined as follows: where H max is the maximum height of the all GSM grid cells inside the ECBilt cell, H sx and H sy are the average peak heights obtained in each row and column of nested cells, and ω 1 is the predefined weight.The silhouette height is then used to calculate the ECBilt cell altitude using: Different combinations of weighting factors ω 1 and ω 2 will draw the gridded terrain analysis toward preserving the peaks (ω 2 = 1) or preserving the valleys (ω 2 = 0), which allows a greater degree of freedom to determine the model terrain analysis.

Ice Mask
Another important consideration in the ice sheet -atmospheric coupling is the variation in ice extent, which changes the albedo calculated by ECBilt and, hence, affects the temperature field over the region and globally.We used the ice thickness field generated by the GSM to create the ice mask needed by ECBilt.To do so, the high resolution ice thickness field is first regridded to the ECBilt coarse resolution grid by using one of the methods mentioned above.Any cell in the resulting grid with more than 30% ice coverage is then assumed to be ice covered.

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 Ocean is sensitive to the freshwater budget at the sites of formation of North Atlantic Deep Water (Rahmstorf, 1995).The numbered circles represent upstream ocean temperature profile sites for ocean-ice coupling.
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. 5b).

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 North Atlantic.This less-saline water affects the upper ocean stratification and thus the strength of deep ocean convection and the AMOC, which in return, 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 crossstrait 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 meter depth) and a closed strait (0.0, 0 meter 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 (given its viscoelastic bedrock response and first order Geoidal correction) has a potentially important role here.

Ensemble Parameters and Sensitivity Analysis
The parameters we chose for our ensemble either represent a crucial physical aspect of the ice sheet -climate interaction (e.g., albedo), or impact the coupling procedure between the ice sheet and the climate (e.g., upscaling method) (Table 2).For each of the ensemble parameters, we applied a Latin-Hypercube scheme to non-uniformly distribute parameter values over the specified interval.The distributions were chosen to increase parametric value density near default parameter values and ensure coverage of a plausible parameter range.
As the context of the model development is glacial inception and deglaciation, we are interested in the ensemble performance for climate metrics which control the growth and decay of the Northern Hemisphere ice sheets during these two stages.
Therefore, we analyze the sensitivity of the coupled model via a focus on summer surface temperature and winter precipitation over land as critical controls of ice sheets growth and decay.To enable comparison against observations, our sensitivity analysis is applied to an ensemble comprised of transient runs covering the historical interval (up to 1980 CE).
As glacial inception and deglaciation trigger at different latitudes in NA and EA, we have divided each continent into diagnostic north and south zones (called "NorthNA", "SouthNA", "NorthEA", and "SouthEA").The sensitivity of the coupled model is tracked for each individual zone.The "NorthNA" and "SouthNA" zones cover latitude ranges of 65-75 • N and 40-60 • N over North America, respectively."NorthEA" and "SouthEA" are defined over 70-80 • N and 55-70 • N latitude bands, respectively.The regional boundaries are illustrated in Fig. 6.
The third column in Table 2 shows the sensitivity of surface temperature and precipitation in the coupled model to changes in each parameter through its range for four different latitudinal bands over NA and EA averaged over the 1950-1980 interval.
For easier comparison, all figures use the same temperature and precipitation scales.The sensitivity to each parameter for the four regions is different for temperature and precipitation.For instance, switching between PD radiative cloud forcing and cloud parameterization strongly affects both temperature and precipitation over all regions, while changing the snow albedo has its strongest impact on EA temperatures and precipitation (Table 2).Each of the ensemble parameters has an impact of at least 4 • C on temperature and/or 1 cm/month on precipitation over the given parameter ranges.We take this as justification for their continued use as ensemble parameters. 5 In the following subsections, we further describe the parameters used in the ensemble simulation.Later, we will show the chosen set of ensemble parameters is adequate for bracketing the relevant (temperature and precip) fields of the climate system.

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 into account the state of the surface (frozen or melting) and the thickness of the snow and ice covers (Goosse et al., 2010).We include all types of snow and ice albedo (ie.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 a cooling effect over all regions with a smaller influence on precipitation, an increase in melting ice albedo causes all regions to get warmer between 4 and 6 • C and increases the winter precipitation.

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 CO2 forcings ranging from 3000 to 5000 years but without 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 spinup 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).

Upscaling
The three different upscaling methods described in section 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 plotted 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.

Precipitation threshold
ECBilt accounts for humidity, and thus precipitable water, only between the surface and the 500 hPa layer.Above 500 hPa, the atmosphere is assumed to be dry, so 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 (H.Goosse et al 2010).This parameter has the largest relative impact for NorthEA (1 cm/month over the parameter range equivalent to 60% of mean).

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-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.The 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 data-set (Rossow, 1996).The total downward and upward shortwave radiative fluxes depend on the transmissivity of the atmosphere, which also relies on the prescribed cloud cover (Goosse et al., 2010).
Given the importance of cloud radiative feedbacks on ice sheet evolution, using a prescribed PD cloud cover for paleoclimate modelling can preclude realistic climate simulations.Therefore, we have added a simple cloud parameterization scheme similar to the precipitation parameterization scheme as described in section 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. 7) through reduced humidity during glacial conditions reducing the cloud cover.

Present-Day Ensemble Results
The fast runtime of the coupled model permitted an initial ensemble of 2000 PD simulations using the fully-coupled GSM-

10
LOVECLIM and varying the model parameters described above.We chose the PD interval to permit comparison of the coupled model output against observational data and to select a better fit sub-ensemble for transient paleo runs.All simulations are spun up using transient forcings (orbital and GHG) for 3000 to 5000 years without the GSM coupled, followed by a transient coupled run ending at year 1980.From here on, we focus on the 500 member sub-ensemble results.Figure 9 shows the Greenland region ensemble mean thickness and standard deviation at PD.The largest ice thickness changes occur at the southern margins of the Greenland ice sheet.Eastern NA ice expansion is concentrated in the high Arctic (Ellesmere Island and adjacent) where PD ice caps exist.

Surface Temperature and Precipitation
The ensemble distribution of the annual mean global surface temperature anomaly with respect to observations shows that the majority of the ensemble members fall within ± 3 • C from observation (grey bars in Fig. 10).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).Green bars in Fig. 10 show the temperature anomaly distribution for the selected 500 members.
Our four latitudinal bands defined in Section 4 (Table 3) provide more relevant temperature metrics for Northern Hemispheric ice sheet contexts.All four regions have higher mean seasonal surface temperatures 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 20

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.They are therefore 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 with observations in winter and summer over NA and EA (Fig. 11a).The ensemble shows good agreement in capturing the maximum zonal velocity, but there is a 10 to 15 degree 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 trough and ridges from the Pacific Ocean to the Atlantic Ocean largely match the reanalysis output within ensemble range (red lines in Fig. 11b) 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.

sea ice
High latitude sea ice acts as an insulator for both heat and moisture between the atmosphere and ocean, the two controlling factors for terrestrial ice sheet surface mass-balance.We use the area and minimum latitude extent of the Northern Hemisphere

Atlantic Meridional Overturning Circulation
The Atlantic Meridional Overturning Circulation (AMOC) transports large amounts of heat and freshwater 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 by at least two AMOC features of the sub-ensemble.The PD ensemble mean AMOC strength is weaker than the reanalysis data from European Center 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. 13a).The ensemble mean shows a slight increase in AMOC strength from 1965 to 1980 CE, which is not seen in ORA-S3.As well, 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. 13a).
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. 13b).

6 Conclusions
We have coupled an Earth System Model of Intermediate Complexity (LOVECLIM) with a 3D thermomechanical coupled ice sheet systems model (GSM).The coupling efficiently captures the main feedbacks between the ice sheet and the atmosphere and ocean models.Our coupled model includes a parameterized sub-shelf melt using marginal ocean vertical temperature profiles, a simple cloud parameterization scheme to improve the radiative forcing representation in the atmosphere model, 10 a dynamical vertical temperature gradient, and a dynamic meltwater runoff routing.We also introduce a new precipitation downscaling scheme that accounts for the change in surface slopes between the coarse resolution climate model grid and the higher resolution ice sheet model grid.Each of the above features has significant impact on modelled ice thickness (shown directly or via changes in temperature and/or precipitation).We have presented a set of ensemble parameters to generate an ensemble of runs that "bracket reality" and have shown that each ensemble parameter has a significant impact on modelled PD regional temperatures and/or precipitation.The new coupled model was subject to a Latin Hypercube parameter sweep of 2000 ensemble simulations for PD boundary and initial conditions.We extracted a sub-ensemble of 500 model runs according to modelled PD Northern Hemispheric ice volume changes.The mean of the sub-ensemble is warm and wet-biased for the Northern Hemispheric ice sheet region.However, the model ensemble still brackets reanalysis precipitation and temperature fields within two (ensemble) standard deviations for all regions and within one standard deviation for half of the regions for the case of temperature.
The ensemble's performance at capturing the seasonal cycle is much better.The ensemble mean difference between summer and winter for all four regions is well within one standard deviation of reanalysis values for both temperature and precipitation (and within one degree C for three of the four regions).This provides some confidence that the model responds adequately to orbital forcing (at least for components that operate on sub-annual time-scales).
The "reality bracketing" criterion is not met for certain features of atmospheric circulation (especially winter-time Rossby wave number) and AMOC strength and variance.Another key limitation of LOVECLIM is the inability to change bathymetry and landmask (aside from the parametrized Bering Strait throughflow).The paleo climate and ice sheet modelling communities would be well served by a modern successor to LOVECLIM for large ensemble glacial cycle time scale contexts.
tines/scripts and modified version of LOVECLIM are freely available upon request.The full LOVECLIM/GSM coupled model is available for collaboration with the authors.
Competing interests.The authors have no competing interests

Table 1 .
Feedbacks included in previous studies between the ice sheet model and the rest of the climate system, compared to the Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1 .
Figure 1.Components of the climate system and interactions between them included in the coupled model.

3. 1 -
Atmosphere to ice At the end of each coupling time-step, the coupler receives climate fields from ECBilt and converts them to monthly-mean values.These fields include: Latitudinal and longitudinal components of wind 6 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.

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

Figure 3 .
Figure 3. Impact of advective precipitation downscaling inclusion in the coupled model; NA ice thickness difference at 110 ka between simulations with and without advective precipitation method.
Figure 5aprovides 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

Figure 4 .
Figure 4. LOVECLIM freshwater discharge mask in a. NA and Gr, and b.EA.Different colors separate regions with same freshwater rates.

Figure 5 .
Figure 5. Impact of meltwater runoff inclusion in the coupled model; a.Total ice volume evolution at glacial inception with (Blue) and without (Red) dynamic meltwater routing, and b.NA ice thickness difference at 110 ka with and without runoff routing.

Figure 6 .
Figure 6.Selected north and south zones over North America and Eurasia for PD sensitivity analysis.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.Table 2. Ensemble 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 mean summer surface temperature and winter precipitation over four selected regions when each parameter is varied independently from its minimum to its maximum value.Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 7 .
Figure 7.Total ice volume in SLE (sea level equivalent) at the last glacial inception with the cloud parameterization (green) and with the PD cloud cover forcing (Red).
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.Given a priority to "bracket reality" and limitations of the component models, we chose to not use climate characteristics for the sub-ensemble filter.Our focus on coupled ice and climate and our choice to avoid bias corrections led to a trial criteria based on ice volume changes (between 1700 and 1980 CE).We first considered a less than 0.1 m SLE (sea level equivalent) change in ice volume requirement for each of the 3 Northern ice sheet regions.But this already fell below our target size of 500 simulations with NA being the problematic region (Fig.8).So we changed the criterion to be the 500 runs with the least amount of ice volume change over each ice sheet.The sub-ensemble NA simulations have ice volume less than 0.15 m SLE and ice volume changes for the other two ice sheets well below 0.1 m SLE (Fig.8).Crucial to our "reality bracketing", there are about 80 sub-ensemble members with ice loss over the given time interval.

Figure 8 .
Figure 8.The distribution of ice volume change over NA, EA, and Gr, between 1700 and 1980 CE in 2000 ensemble runs.Green bars represent the selected 500 ensemble members, and red bars represent the rest.

Figure 9 .
Figure 9. Greenland ice thickness ensemble a. mean, and b. standard deviation at PD.

Figure 10 .
Figure 10.The distribution of global annual mean surface 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.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-277Manuscript under review for journal Geosci.Model Dev. Discussion started: 2 January 2018 c Author(s) 2018.CC BY 4.0 License.regions have a mean difference between summer and winter ensemble temperatures within a degree C of that of the reanalysis climatology.

Figure 11
Figure 11.a. Zonally average of the zonal component of the 200 hPa wind velocity, and b. meridionally average of the meridional component of the 200 hPa wind velocity over the Northern Hemisphere.Filled areas show the model ensemble mean and the two standard deviation range.Dashed lines represent observational data.Blue is for winter and red for summer.

Figure 12 .
Figure 12.Maximum (March: blue) and minimum (September: red) sea ice area ensemble mean ± one standard deviation.The vertical lines represent observational 1981-2010 March and September mean sea ice area within one standard deviation(Walsh et al., 2015).

Figure 13
Figure 13.a. Maximum AMOC strength at 26 • N between 1966 and 1980 CE. Black solid line: ensemble mean; dark blue area: ensemble mean ± one standard deviation; light blue area is bounded by the simulations with the maximum and minimum AMOC strength (time averaged); dashed line: ORA-S3 ((Balmaseda et al., 2008)).b, AMOC stream-function mean (filled colors) and ensemble standard deviation (contour lines).

Table 3 .
Summer surface temperature and winter precipitation averaged over four latitudinal bands for summers and winters of