Effect of orographic gravity wave drag on Northern Hemisphere climate in transient simulations of the last deglaciation

Long transient simulations of the last deglaciation are increasingly being performed to identify the drivers of multiple rapid Earth system changes that occurred in this period. Such simulations frequently prescribe temporal variations in ice sheet properties, which can play an important role in controlling atmospheric and surface climate. To preserve a model’s standard performance in simulating climate, it is common to apply time dependent orographic variations, including parameterised sub-grid scale orographic variances, as anomalies from the pre-industrial state. This study investigates the causes of two abrupt climate change events in the Northern Hemisphere extratropics occurring between 16 and 14 thousand years ago in transient simulations of the last deglaciation from the Hadley Centre coupled general circulation model (HadCM3). One event is characterized by regional Northern Hemisphere changes comprising a centennial scale cooling of ~ 10 °C across Fennoscandia followed by rapid warming in less than 50 years as well as synchronous shifts in the Northern Annular Mode. The second event has comparable but temporally reversed characteristics. Sensitivity experiments reveal the climate anomalies are exclusively caused by artificially large values of orographic gravity wave drag, resulting from the combined use of the orographic anomaly method along with a unique inclusion of transient orography that linearly interpolates between timesteps in the ice sheet reconstruction. Palaeoclimate modelling groups should therefore carefully check the effects of their sub-grid scale orographic terms in transient palaeoclimate simulations with prescribed topographic evolution.


Introduction
Long transient simulations of palaeoclimate are increasingly being performed to better understand the drivers of the sequence of climate events that occurred during the last transition from cold glacial to warm interglacial conditions (e.g. Liu et al. 2009;Menviel et al. 2011;Roche et al. 2011;Bethke et al. 2012;Gregoire et al. 2012;Obase and Abe-Ouchi 2019). Known as the last deglaciation, this transition was initiated by a gradual increase in boreal summer insolation ~ 23 thousand years before present (ka BP; Berger 1978), causing the vast ice sheets spanning much of the Northern Hemisphere continents (Dyke et al. 2002;Lambeck et al. 2014;Hughes et al. 2016) to begin melting (Gregoire et al. 2011). The subsequent increase in atmospheric greenhouse gases (e.g. Loulergue et al. 2008;Schilt et al. 2010;Bereiter et al. 2015) then became the dominant driver of much of the remaining deglaciation (Gregoire et al. 2011).
A particularly compelling characteristic of the last deglaciation is a series of rapid earth system changes, including abrupt (decadal to centennial) hemispheric-scale warmings and coolings of more than 5 °C (de Beaulieu and Reille 1992;Severinghaus and Brook 1999;Lea 2003;Buizert et al. 2014), sudden reorganisations of basin-wide ocean circulation (e.g. Roberts et al. 2010;Ng et al. 2018), and jumps in sea level of tens of metres (e.g. Deschamps et al. 2012;Lambeck et al. 2014). As such, the last deglaciation presents a particularly useful case study for examining and testing our models' abilities to simulate the ice-ocean-atmosphere interactions that operate under warming climates with melting ice sheets and abrupt (human timescale) changes. Yet, despite decades of research into this most recent period of Earth's geological past, the precise chain of events that followed and their underlying causal mechanisms remain enigmatic (e.g. see discussions in Bethke et al. 2012;Ivanovic et al. 2018;Obase and Abe-Ouchi 2019). Thus, under the auspices of the Palaeoclimate Modelling Intercomparison Project Phase 4 (PMIP4), a protocol was published providing information on how to run transient simulations of the last deglaciation  to facilitate a coordinated effort for cross-model comparison tapping into the full range of available climate and earth system models as a resource for examining uncertainties in deglacial forcings, trigger mechanisms, and dynamic feedbacks.
One notably complex aspect of simulations of the last deglaciation is the common requirement to prescribe the global distribution of continental ice sheets, including timeevolving orography, as ice mass loading changes. Changes in ice sheet elevation are already considered to be the most important factor to influence atmospheric circulation and sea level pressure (Pausata et al. 2009(Pausata et al. , 2011Löfverström and Lora 2017;Gregoire et al. 2018), both potentially deriving from the sub-grid scale parameterisations of anisotropic variations in orographic height (known as orographic variance) within model simulations. For example, Gong et al. (2015) and Gregoire et al. (2018) show through model simulations that higher elevations of the Laurentide and Greenland ice sheets correspond to strong surface wind stress curl over the North Atlantic. Some more recent earth system models demonstrate similar behaviour with the inclusion of dynamic ice sheet models (e.g. see Lofverstrom and Liakka 2018;Lofverstrom et al. 2020;Ackermann et al. 2020;Gregory et al. 2020), and some groups have developed dynamiccontinent components within this framework (Meccia and Mikolajewicz 2018). However, ice sheets and their associated geographical features are still mostly externally prescribed in palaeoclimate simulations. If a prescribed ice sheet is required, the PMIP4 protocol for simulating the last deglaciation provides a choice of two global ice sheet reconstructions spanning from 26 ka BP to present, ICE-6G_C VM5a (henceforth 'ICE-6G_C'; Argus et al. 2014;Peltier et al. 2015) and GLAC-1D (Tarasov and Richard Peltier 2002;Tarasov et al. 2012;Briggs et al. 2014;Ivanovic et al. 2016).
As well as orography directly impacting the atmospheric circulation, it can also have an indirect effect via parameterisations of sub-grid scale orographic gravity waves. These waves are generated by flow over orography and transport momentum from their source to where they dissipate or are absorbed, thereby imparting drag on the mean flow (McFarlane 1987). It is necessary to parameterise this momentum transport in atmospheric models due to the spatial scale of the waves being much smaller than the model grid (Kim et al. 2003) and because omitting the drag leads to climatological biases in circulation (Shaw and Shepherd 2007;Sandu et al. 2016;Pithan et al. 2016), such as excessively strong westerlies in the midlatitude Northern Hemisphere (Kim and Arakawa 1995). Since the first implementation of an orographic gravity wave drag parameterisation in a large-scale model by Boer et al. (1984), multiple other schemes have been developed, for example by Chouinard et al. (1986), Palmer et al. (1986), McFarlane (1987, and Lott and Miller (1997). Most modern general circulation and Earth system models still use variations of these early schemes; for example, the Met Office Unified Model uses a development of the Palmer et al. (1986) scheme but with some improvements to better represent low and mid-level drag (Gregory et al. 1998), and the Lott and Miller (1997) scheme is referenced for the CMIP6 Last Glacial Maximum experiment as described by Kageyama et al. (2017). Orographic variance also impacts form drag, increasing surface friction over rough orographic features. In the Hadley Centre Coupled Model version 3 (HadCM3; Valdes et al. 2017), this is parameterised by increased effective roughness length (Gregory et al. 1998) and depends on the sub-grid scale orographic variance. Despite the abundance of evidence for the importance of orographic gravity wave drag for simulating the general circulation, few studies have examined the interaction between palaeoclimate orographic variations and parameterized orographic gravity wave drag. This stands in sharp contrast to other aspects of palaeoclimate surface and boundary conditions such as ice sheet/shelf extents, coastlines (including inland seas), river routing and bathymetry (e.g. Ivanovic et al. 2016;Kageyama et al. 2017;Menviel et al. 2019).
In this study, we present transient simulations of the last deglaciation using HadCM3 (Valdes et al. 2017), following the PMIP4 protocol  with the use of ICE-6G_C to define the palaeo ice sheets (Peltier et al. 2015). Unique to these simulations, orography is prescribed as a transient term derived from linear interpolation between the time steps in the reconstruction (every 1000 years 23-21 ka BP and every 500 years thereafter). We show that very small inconsistencies in the calculation of orographic sub-grid scale variances can produce intervals of artificially high orographic gravity wave stress over the North American ice sheet, with profound impacts on the Northern Hemisphere extratropical circulation and regional climates. This phenomenon reaches an extremity at 15.5 ka BP and 14.5 ka BP, with both cases traceable to a spurious behaviour of the orographic gravity wave drag scheme for the prescribed palaeo-orographic boundary conditions.
In the following text, we outline the model set-up and sensitivity experiments used to isolate the cause of this event (Sect. 2), present and discuss our results (Sect. 3), and conclude with a cautionary note that in future simulations with transient palaeogeographies, special care must be taken in the model set-up of sub-grid scale parameterisations. For brevity, we focus the presented analysis on the 15.5 ka BP event.

Model description
We performed transient simulations of the last deglaciation with the HadCM3 coupled ocean-atmosphere-vegetation general circulation model (GCM; Gordon et al. 2000;Pope et al. 2000) with modest modifications described by Valdes et al. (2017). This model, often used to simulate palaeoclimates, has an atmospheric resolution of 2.5° latitude by 3.75° longitude with 19 vertical levels starting at the surface and ending at 10 hPa. The ocean's horizontal resolution is 1.25° by 1.25° with 20 vertical layers from the surface of the ocean to ~ 5500 m depth with maximum resolution at the surface. Coupled to the atmospheric GCM, vegetation is represented by the dynamic vegetation model Top-down Representation of Interactive Foliage and Flora Including Dynamics (or TRIFFID) and linked to the land surface with the Met Office Surface Exchange Scheme or MOSES 2.1. Despite its lower resolution, the performance of HadCM3 in simulating mean climate is comparable to other more recent models (Flato et al. 2013;Reichler and Kim 2008), and importantly, it is sufficiently computationally efficient to make it feasible to perform long integrations, such as the 25,000-year long simulations analysed here, within several months.

Orographic gravity wave drag parameterisation in HadCM3
The parameterisation of orographic gravity wave drag in HadCM3 is based upon the scheme implemented in the Met Office Unified Model. The theoretical basis for the gravity wave drag scheme comes from Shutts (1990) and replaces the previous scheme based on Palmer et al. (1986). The four orographic variance terms that relate to the gravity-wave stress are the standard deviation of the sub-grid scale orography ( ), the stress component on the x-plane in the x-direction ( xx ), stress component on the x-plane in the y-direction ( xy ), and the stress component on the y-plane in the y-direction ( yy ). As is customary in palaeoclimate modelling (see discussion in Sect. 1), they are calculated using the data provided by ICE-6G_C as anomalies from pre-industrial, rather than as absolute fields. These four components contribute to the hydrostatic surface stress as given by where, s is the surface values of the density, N s is the Brunt-Vaisala frequency, U 0 is the component of the wind in the direction of the surface stress, ̂ is a typical wavelength for the spectrum of gravity-waves being parameterised, is the direction if the surface layer wind relative to westerly, is the standard deviation of the sub-grid scale orography, xx , xy , and yy are the squared gradients of the elevation of the sub-grid scale orography.
The values of ĥ 1 and ĥ 2 depend on the low-level Froude number and are calculated as: where, and are constants to be specified, determining the regime for a given Froude number (Eq. (5), below), N s is the Brunt-Vaisala frequency, U is the measure of the wind speed of the low-level flow (Walters et al. 2019). The Froude number is defined bŷ is defined by Shutts (1990) as: where, k u is the highest wavenumber lying outside the range of trapped lee waves and k 1 is the lowest unresolved wavenumber of the model.

Last deglaciation
Simulations of the last deglaciation were run with boundary conditions following the PMIP4 forcing protocol described by Ivanovic et al. (2016). These simulations of the entirety of the last deglaciation began from the end of a previous equilibrium simulation for 23-ka (as described and utilised by Morris et al. 2018) and continued with transient boundary conditions/climate forcings until 2000 years into the future (3950 CE). For the set of simulations presented here, we use the ICE-6G_C ice sheet reconstruction (Argus et al. 2014;Peltier et al. 2015), which provides the input data for prescribing ice sheet topography and extent, as well as associated boundary conditions (orography, bathymetry, land-sea mask, ocean-mask, and global mean ocean salinity) at 1000 years intervals from 23 to 21 ka BP and 500-year intervals from 21 to 0 ka BP. Orbital parameters (as per Berger 1978)  Ice sheet meltwater inputs to the ocean follow the 'meltuniform' scenario proposed by Ivanovic et al. (2016), where we prescribe global mean ocean salinity according to a linear interpolation of the change in ice volume between the ICE-6G_C time-steps. Any difference between this target global mean salinity (from ICE-6G_C) and simulated global mean salinity is redistributed evenly across the entire volume of the ocean at each timestep. There remains significant uncertainty in the history of freshwater forcing from deglaciating ice sheets (e.g., discussion by Ivanovic et al. 2016), and we ran several other meltwater forcing and meteoric river routing scenarios to test its impact on our simulations, including versions of the alternative 'melt-routed' scenario proposed by Ivanovic et al. (2016). In terms of the processes presented in this study, all of these simulations yielded the same results, thus for succinctness and clarity, we present the analysis of just one 'melt-uniform' freshwater-forcing scenario, henceforth referred to as degl. Every time we make an adjustment to the model's boundary conditions (i.e. when implementing changes in land-sea mask, bathymetry, and ocean-mask at the 1 ka intervals 23-21 ka, and 500-year intervals thereafter), we branch off the simulation to produce an overlapping 50-year segment where these boundary conditions are not changed. Thus, by comparing the overlapping 50 years, we can verify the effects of introducing the new paleogeographic configurations, and evaluate the magnitude of any resulting climateshock (e.g. as discussed/shown in the climate forcing used by Gregoire et al. 2012Gregoire et al. , 2016.

Orography boundary condition
To incorporate the changes in ice sheet elevation into the climate model, it is common to apply these reconstructions as anomalies from the pre-industrial state. This is advocated for most of the simulations supported by PMIP4. In such cases, the orographic height for a time period is equal to the present-day orography (as in the "control" model simulation) plus the change of height as defined by the appropriate protocol. Thus, the total elevation at any given time period is the pre-industrial elevation plus the change in elevation given by the reconstruction (suitably regridded onto the model grid). This methodology preserves the model's present-day configuration and hence does not require rerunning the pre-industrial control.
We also must calculate several components of the subgrid scale orographic variations, and this is more problematic. An ice sheet will typically be relatively smooth in its interior but with large sub-grid scale variability at its edges. Hence in general, we would expect the variances to be smaller than present day in the interior of an ice sheet but large around the edge. It is further complicated in high mountain regions where the mountain peaks may extend beyond the ice sheet itself. The fundamental problem is that high spatial resolution orographic data (typically 10 km or less) is required to calculate sub-grid orographic variances and these are not available for paleo time periods. Some paleo boundary conditions may be presented at sufficiently high resolution (e.g. ICE-6G_C is 10 min resolution and GLAC-1D is 0.5°) but in practice there is often relatively little variance at these fine resolutions. Hence it is not possible to accurately calculate the variances from any existing paleo datasets, and modelling protocols vary in terms of prescriptive. The PMIP4 LGM protocol (Kageyama et al. 2017) recommends applying the 10-min orographic differences to an existing presentday orography field, and subsequently smoothing the data using a Gaussian filter, but this will potentially result in an underestimate of total variability. Alternatively, the PMIP4 last deglaciation protocol ) specifies using the anomaly method. Since we are performing PMIP4 last deglaciation simulations, we choose to follow this latter approach, and thus determine all variance terms by calculating the change in orographic variance and adding it to the existing pre-industrial values. This will also likely underestimate change because if the anomaly in 2 , xx , or yy results in the totals becoming less than zero, they are reset to zero, but it otherwise preserves the total variability.
The total orography elevation and variances (preindustrial + anomaly) are then applied to the model. The input data is at every 500 years. Rather than having abrupt elevation (and variances) changes every 500 years, every month we linearly time interpolate all orographic fields so that all vary gradually and smoothly between the 500-year tie-points. An example of the evolution of these fields are given in Fig. 5b, which shows the time evolution of in 2 and the maximum of xx , or yy for a specific locality. Note that the variations of the max( xx , yy ) is linear, whereas the 2 varies quadratically because the time interpolation is linear in σ.

Sensitivity experiments
In addition to the full-deglaciation simulations described above, we ran multiple shorter sensitivity experiments to isolate the cause of the abrupt climate change simulated at ~ 15.5 ka (Table 1). These simulations were performed with the atmospheric GCM only to remove the potential influence of the ocean. Therefore, for the ocean boundary conditions, the sensitivity experiments use the equilibrium climate produced by a steady-state simulation of 16 ka (equivalent to the one used by Morris et al. 2018). Our control simulation for these tests, known as 'a-control', is set with 15.5 ka BP initial boundary conditions and integrated for 100 years. We ran an additional control experiment, known as '15.45_control', which was set with the same initial climate conditions as a-control and integrated for 100 years, but has 15.45 ka BP boundary conditions imposed. This experiment provides a baseline for assessing the change in the surface climate after the abrupt 15.5 ka BP warming event and further proves that the ocean does not have a major influence on the anomalies. 15.45_control forms the basis for nineteen additional simulations that systematically vary the atmospheric boundary conditions with respect to the control, homing in on sub-grid scale The longitudinal cut-off point for delineating western and eastern North America (for East_NA_orog_constant and West_NA_orog_constant) is defined at 105° W. For the four orographic variance terms, σ is the standard deviation of the sub-grid scale orography (sd_orog_only), σ xx is the stress component on the x-plane in the x-direction (GWxx_orog_only) σ xy is the stress component on the x-plane in the y-direction (GWxy_orog_ only), σ yy is the stress component on the y-plane in the y-direction (GWyy_orog_only); see Sect. orographic variance (defined in Sect. 2.1.1) in increasing detail. All nineteen of these supplementary sensitivity simulations utilise the same initial conditions as a-control and are integrated for 100 years, see Table 1 for details.

Calculation of Northern Annular Mode index
We examine the Northern Annular Mode (NAM) which is the dominant mode of variability in the Northern Hemisphere extratropical circulation (Thompson and Wallace 2001). Changes in the NAM represent fluctuations in the distribution of atmospheric mass between mid-and high-latitudes, with a more positive [negative] NAM index indicating a stronger meridional atmospheric pressure gradient from zonally lower [higher] pressure in north to higher [lower] pressure in the south. Following Li and Wang (2003), we calculate the NAM index as the difference in zonal mean sea level pressure at the model grid boxes closest to 35° N and 65° N. Statistical significance of the sensitivity runs was confirmed by performing two-sided t-tests at the 95% confidence level.

Results and discussion
In the degl simulation, an unexpected abrupt climate change occurs across the Northern Hemisphere between 16 and 15.45 ka BP that does not correspond directly to changes in solar orbit, atmospheric greenhouse gas concentrations, ice extent, bathymetry/coastlines or freshwater forcing from ice sheet melting. The changes in climate over this period are particularly prominent in the North Atlantic and European region in boreal winter, with surface air temperatures over Fennoscandia and the Nordic Seas regions cooling by ~ 10 °C within around 100 years and then suddenly increasing by a similar magnitude in less than a decade (Fig. 1c). The December-January-February mean surface temperature anomalies over Greenland and the Labrador Sea also show significant changes of ~ 2-3 °C in opposing signs. These changes can be compared to the small (0.25 °C) decrease in global mean surface air temperature in boreal winter between 16 and 15.45 ka BP, indicating the strong locality of the climate patterns. The surface temperature changes are coincident with significant increases in precipitation (Fig. S1) and marked variations in the NAM index (Fig. 1d). During the period of Fennoscandian cooling, the winter NAM trends towards a more negative state (i.e., weakened meridional pressure gradient), decreasing by around 5 hPa, and then increasing very rapidly over a period of around a decade back to a state equivalent to that around 15.8 ka BP. Similar behaviour of a smaller magnitude is evident in boreal autumn and spring. The change in winter NAM index over this period is substantially larger than the magnitude of inter-decadal variability, at ~ 4σ above the decadal average, and is differentiated from the NAM variability over the remainder of the simulation (though similar variations are seen around 14.5 ka BP when another period of rapid regional climate change is simulated-see Sect. 1).
The spatial patterns of surface climate anomalies during the transition in climate state around 15.5 ka BP are shown in Fig. 2. During the period of Fennoscandian cooling, there is a weakened meridional pressure gradient across the northern extratropics at most longitudes (Fig. 2d). This shows a dominant centre of action in the North Atlantic, resembling the North Atlantic Oscillation (NAO), but interestingly, the dipole in mean sea level pressure (MSLP) also extends across the Eurasian continent, which is not a canonical feature of the NAM (Thompson and Wallace 2001). This pattern represents a redistribution of atmospheric mass from midlatitudes to high latitudes. During the period of rapid Fennoscandian warming, the pattern of MSLP anomalies is very similar, but opposite in sign (Fig. 2e). Figure 2a and b shows that in the period of a negative NAM trend, there is cooling over western Europe and Fennoscandia, warming over eastern Canada and the USA, and warming over North Africa. This is the canonical pattern of surface temperature anomalies associated with the NAO through modified thermal advection (Thompson and Wallace 2001). Across the Eurasian continent, there is cooling at high latitudes and warming at low latitudes, which reflects reduced poleward heat advection (Fig. 2a). The regional temperature anomalies Hatching denotes regions where the anomalies are not statistically significantly different from a-control at the 95% confidence level using a Student t-test are reversed in sign in the later period with a rapid increase in NAM (Fig. 2b; although the temperature anomaly over the Western Atlantic is not statistically significant within the natural variability, and therefore is not shown). Hence, the regional surface temperature anomalies appear to be a consequence of the anomalous atmospheric circulation. This is consistent with the radiative forcing being small in this interval (Fig. 1e), and with the changes in global mean surface temperature being small.
The timing of the abrupt changes in climate (Figs. 1, 2) provides a strong indication that this event is not directly driven by ice sheet volume change since there are no coincident large changes in terrestrial ice volume (furthermore, note that there is no surface freshwater forcing from ice sheet meltwater in the PMIP4 'melt-uniform' scenario employed here). Contrary to many of the previous mechanisms identified for driving abrupt transitions in climate in deglacial/ glacial simulations, where much larger changes occur, (e.g. Knutti et al. 2004;Liu et al. 2009;Li et al. 2010;Dokken et al. 2013;Obase and Abe-Ouchi 2019), the differences in sea ice concentration, extent and thickness, and the Atlantic meridional overturning circulation over the period of interest (as demonstrated in Fig. S2) are negligible in comparison to the magnitude of the climate change, suggesting that the underlying processes affecting the atmospheric circulation do not originate in the ocean. This is confirmed by the atmosphere-only sensitivity simulations spanning 15.5 ka and 15.45 ka (see Table 1 and Sect. 2.2.3; a-control and 15.45_control), which produce the same anomalies in surface climate despite the absence of a coupled ocean. We next seek to explain the source of the anomalous northern extratropical circulation in the degl simulation.
To isolate the causal mechanism, we performed and analysed the 19 additional sensitivity simulations described in Table 1, which hold certain boundary conditions fixed at 15.5 ka BP values, while incrementally updating others to match conditions at 15.45 ka BP. Results from these sensitivity experiments were compared to a-control, which uses the full suite of 15.5 ka BP boundary conditions.
The sensitivity experiments that test the effect of radiatively active gases and orbit (RAGS_constant), land surface properties (TRIFFID_constant), and sea surface temperature and sea ice (SST_ice_constant) all show the same North Atlantic climate anomalies as in degl (as demonstrated by Figs. 2, 3a). Since orography is the only component not reverted to 15.5 ka BP conditions in these simulations, we determine that the definition of orography must be driving the climate anomalies in this period.
Further sensitivity experiments were performed, iterating from SST_ice_constant (where only orography is updated to 15.45 ka BP), to determine which component(s) of the orography fields drive the event. Individual orographic definitions and locations were reverted to 15.5 ka BP conditions, as in a-control, in simulations orog_var_constant, all_but_US_Asia_orog_constant, Fenn_orog_constant, UK_Iceland_orog_constant, Greenland_orog_constant, East_NA_orog_constant, and West_NA_orog_constant. Although small differences occurred from the changes in orography in each of the experiments, the main features of the abrupt climate event (e.g. Fig. 2b) were present in all simulations except for orog_var_constant (Fig. 3b), when all orographic variances are held at 15.5 ka conditions; this demonstrates that the orographic variances are driving the simulated changes of main interest.
To determine precisely which aspect of the orographic variances is key for driving the anomalous atmospheric circulation, further experiments were performed: SD_orog_ only, GWxx_orog_only, GWxy_orog_only, GWyy_orog_ only, silhouette_orog_only, and peak-to-trough_orog_only (Table 1; as displayed in Fig. 4). In each simulation, all cally significantly different from a-control at the 95% confidence level using a Student t-test variances were held at 15.5 ka BP conditions, as in a-control, except for the individual parameterisation of orographic variance being tested ( , xx , xy , and yy , silhouette, and peak-to-trough, respectively), which is set to 15.45 ka BP conditions as in 15.45_control. As shown in Fig. 4b-f, the adjustments to xx , xy , and yy , silhouette, and peak-totrough orographic variance all instigate only minor changes in surface temperature that are largely statistically insignificant. Only the standard deviation of the orography ( ) has a strong, and statistically significant, effect on surface air temperature (Fig. 4a), reproducing a similar pattern of climate changes as observed in degl.
Based on these results, we conclude that the parameterisation of sub-grid scale orographic variance has a strong influence on climate during the 16-15.45 ka period as simulated in the HadCM3 model. Further simulations sd_orog_ all_areas, sd_orog_NA_only, and sd_orog_one_grid_pt_ only (Table 1) are used to pinpoint the specific location(s) that influence the climate evolution. Surprisingly, the abrupt climate event as evident in the degl simulations and Fig. 2b is closely reproduced in the sd_orog_one_grid_pt_only simulation. In this final sensitivity simulation, only the standard deviation of the orography of one grid point in northwestern North America (62.5° N, 123.75° W) is changed to 15.45 ka BP conditions, everywhere else the conditions correspond to those in a-control (Fig. 5a). Figure 5b depicts a time series of the orographic variance and hydrostatic stress at the grid box centred over 62.5° N, 123.75° W during the time of the abrupt climate event shown in Figs. 1 and 2. This shows the time evolution of σ 2 , the largest of σ xx and σ yy and the ratio of the two. During the 500-year interval between 16 and 15.5 ka BP, the orographic variance decreases, reaching zero just before 15.5 ka BP. This has a profound impact because the orographic variance provides the denominator of the hydrostatic stress equation (Eqs. 1, 2). Thus, as the orographic variance approaches zero, the hydrostatic stress becomes artificially and unphysically large at a rapidly accelerating rate. This induces a southerly shift in the track of the polar jet stream, which causes the centennialscale transient reduction in the NAM index and associated cooling over Fennoscandia. When the orographic variance reaches zero at 15.5 ka BP, the orographic gravity wave drag parameterisation is configured to respond by turning off at that location, allowing the westerly winds to quickly resume alongside the recovery of the polar jet stream track, resulting in an instant rebound of the NAM, which drives the abrupt temperature change of ~ 10 °C in Fennoscandia.
Even though the prescription of the orographic variance is not seasonal, there is a pronounced intra-annual asymmetry in the resulting climatic signal, as shown in Fig. 1a-d, with the greatest anomalies occurring in the boreal winter and negligible impact in the June-July-August. In the boreal summer, the low-level westerly winds are located farther south (as illustrated in Fig. S3) causing the U 3 term in Eqs. (1) and (2) to be small, decreasing the overall magnitude of the hydrostatic stress and hence reducing the impact of the temporal evolution in the orographic variance centred at 62.5° N, 123.75° W described above.
The described unphysical growth of the hydrostatic stress is a direct consequence of using the commonly employed is silhouette_orog_only, and f is peak-to-trough_orog_only. Hatching denotes regions where the the anomalies are not statistically significantly different from a-control at the 95% confidence level using a Student t-test anomaly-from-present method for deriving palaeo-environmental orographic boundary conditions, alongside the application of linear interpolation to smooth fields between time points. The coincidence in location of the unphysical term and the path of the winter low-level westerlies induces the profound effect on atmospheric climate observed in our simulations (with negligible effect on the ocean). We have verified that the similar climate event simulated at 14.5 ka BP in degl is also explained by the same technical phenomenon, although from a different change in ice sheet geometry. We do not observe this process at any other time in the 25 ka BP simulation.
The timing of the events at 15.5 ka BP and 14.5 ka BP coincides with key points of deglaciation in the expression of the Laurentide-Cordilleran ice sheet saddle collapse (Gregoire et al. 2012) in ICE-6G_C. Between 16 and 14 ka in the ice sheet reconstruction, ~ 1500 m of ice is lost in this region. However, at the grid-point of interest (62.5° N, 123.75° W), the ice sheet thins by ~ 500 m in 500 years between 16 and 15.5 ka BP. Such thinning rates are not unusual in ablation areas over the deglaciation. That such a common change in ice thickness could induce a doubling of the wind stress and such rapid, high amplitude swings in the NAM (producing ~ 10 °C temperature change over Fennoscandia, ~ 3 °C over Greenland and ~ 3 °C over the Labrador Sea, all in a few years, for example) provides a stark but hitherto un-noticed and thus highly useful demonstration of the importance of details in the technical formulation of orographic variance and their relationship to gravity wave drag.
Notably, in the simulations presented here, none of the parameterisations of the sub-grid scale orographic variance terms are out of the ordinary and only smooth transitions between model time steps occur (in part shown by Fig. 5b). Furthermore, the extreme behaviour of the orographic gravity wave drag parameterisation scheme identified at two points in the deglaciation simulation can be hidden at times when other conditions are varying or when extreme climate events are expected/forced to happen (e.g. the Bolling Warming, Meltwater Pulse 1a, Younger Dryas, Dansgaard-Oeschger events more broadly). We determine that the climate transitions seen in our simulations would have been less likely to occur if the orographic fields had been calculated from absolute fields rather than as anomalies from pre-industrial. In addition, because of when they occur in our simulations, we know that they would also not have taken place if orography did not evolve smoothly between each model timestep. These prerequisites, along with the known influence of orographic height on the surrounding atmosphere and surface climate, highlight the importance of carefully checking the evolution and effects of the subgrid scale orographic terms, especially, for example, when attempting to diagnose the results of a large change in the ice sheets, which may otherwise hide the behaviour found in the degl simulations.

Conclusions
This study investigated the drivers of two unexpected abrupt climate changes at ~ 15.5 ka BP and ~ 14.5 ka BP in a long transient simulation of the last deglaciation performed with the HadCM3 model. The climate anomalies comprise large, rapid changes in the Northern Hemisphere extratropical circulation and associated regional patterns Results of the sd_orog_one_grid_pt_only sensitivity experiment. a Demonstrates the change in surface air temperature between the sensitivity run and a-control. The single grid cell with 15.45 ka BP conditions is indicated by the cross (×; 62.5° N and 123.75° W), conditions everywhere else are kept the same as in a-control. Hatching denotes regions where the anomalies are not statistically significantly different from a-control at the 95% confidence level using a Student t-test. b Shows the change in orographic variance in the marked grid cell from 16 to 15.4 ka, where the solid line is the subgrid scale standard deviation (σ 2 ), the dashed line is the maximum value of the squared gradients of the elevation of the sub-grid scale orography (max(σ xx , σ yy )) and the dotted line is the ratio between the two ((max(σ xx , σ yy ))/σ 2 ). All are at 62.5° N and 123.75° W and the squared gradient (σ xx , σ yy ) have been multiplied by 10 5 for plotting of surface temperature change. These events are shown to be a consequence of the orographic gravity wave drag parameterisation scheme in the model, resulting from the uniqueness of the prescription of the orographic boundary conditions in these simulations, and the fact that these technical anomalies coincide with the track of the lower level winter westerlies. The combined use of an anomaly method, where sub-grid scale parameterisations are calculated from the data provided in the ice sheet reconstruction as a difference from the pre-industrial, alongside a transient-evolving orography and more frequent than usual time steps in the reconstruction (every 1000 years 23-21 ka BP and every 500 years thereafter) leads to an unrealistic amplification of the hydrostatic stress in the orographic gravity wave drag parameterisation scheme leading to pronounced impacts on the large-scale atmospheric circulation. This study provides an important outlook on the configuration of palaeoclimate simulations and specifically the parameterisation of sub-grid scale orographic variance and its effect on the simulated surface climate. More thorough communication between modellers utilising current models and palaeoclimate modellers using their earlier counterparts could provide more understanding of how parameterisations are implemented, as often the parameterisation set-ups are the same. We have shown that even very small and seemingly physically realistic changes in the orographic variance can cause a surface climate change of ~ 10 °C of local seasonal warming and cooling, a large response to a technical issue that could, in other cases, present itself more modestly. Because the last deglaciation has become a time period of increasing interest due to the series of rapid earth system changes that are scattered throughout, it is crucial to ensure parameterisations are behaving in the manner intended to develop rational and comprehensive hypotheses. As such, these results highlight the importance of checking a model's sub-grid scale orographic variance terms and their effect on the simulated climate when setting up a new geographical configuration in palaeoclimate simulations with prescribed topographic evolution.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.