High-resolution LGM climate of Europe and the Alpine region using the regional climate model WRF

In this study we present a series of sensitivity experiments conducted for the Last Glacial Maximum (LGM, ∼21 ka) over Europe using the regional climate Weather Research and Forecasting model (WRF). Using a four-step two-way nesting approach, we are able to reach a convection-permitting horizontal resolution over the inner part of the study area, covering central Europe and the Alpine region. The main objective of the paper is to evaluate a model version including a series of new developments better suitable for the simulation of paleo-glacial time slices with respect to the ones employed in former studies. The evaluation of the model is conducted against newly available pollen-based reconstructions of the LGM European climate and takes into account the effect of two main sources of model uncertainty: a different height of continental glaciers at higher latitudes of the Northern Hemisphere and different land cover. Model results are in good agreement with evidence from the proxies, in particular for temperatures. Importantly, the consideration of different ensemble members for characterizing model uncertainty allows for increasing the agreement of the model against the proxy reconstructions that would be obtained when considering a single model realization. The spread of the produced ensemble is relatively small for temperature, besides areas surrounding glaciers in summer. On the other hand, differences between the different ensemble members are very pronounced for precipitation, in particular in winter over areas highly affected by moisture advection from the Atlantic. This highlights the importance of the considered sources of uncertainty for the study of European LGM climate and allows for determining where the results of a regional climate model (RCM) are more likely to be uncertain for the considered case study. Finally, the results are also used to assess the effect of convection-permitting resolutions, at both local and regional scales, under glacial conditions.

1 Introduction

The Last Glacial Maximum (LGM) is known as the time at which the Northern Hemisphere ice-sheet volume reached its maximum during the last glacial period, at around 21 ka (Clark et al., 2009). The radiative forcing of the LGM was considerably different than the preindustrial (PI) period. Atmospheric greenhouse gas (GHG) concentrations were lower at LGM than at PI. At the same time, changes in the orbital configuration of the Earth around the Sun induced a different seasonal pattern of incoming insolation during the LGM, with lower values over the higher latitudes of the Northern Hemisphere during summer, with respect to the present day (Berger, 1978; Loulergue et al., 2008; Kageyama et al., 2017; Bereiter et al., 2015). The direct effect of different radiative forcing led to important changes at the global scale in the climate of the LGM compared to the present day, with overall cooler and drier conditions. Syntheses of proxy data and model outputs indicate that annual global temperatures were in the range of 4 to less than 8 °C lower than their PI values (Osman et al., 2021; Hargreaves et al., 2011). At the same time, feedback mechanisms, such as the ice-albedo effect, land cover changes, and ice-sheet expansion, played an important role in modulating the climate of the LGM, not only at the global scale but also at a regional and local level. For example, a large extension of northern hemispheric ice sheets had a strong impact on the large-scale circulation, with a different downward impact over different parts of Europe (Hofer et al., 2012b, a; Löfverström et al., 2014; Merz et al., 2015; Löfverström and Liakka, 2016; Löfverström et al., 2016; Ludwig et al., 2016; Kageyama et al., 2017).

The study of the climate of the LGM offers the opportunity to better understand processes and feedbacks of the climate system that have no analogue for the present day and the near-future (Kohfeld and Harrison, 2000; Kageyama et al., 2017; Raible et al., 2020). For this reason, the LGM represents a unique opportunity for evaluating the response of climate models to an extreme change in climate forcing, improving model reliability and contributing to model development (Kohfeld and Harrison, 2000). Consequently, the LGM has been one of the main target periods of paleoclimate modeling studies, with the first attempts at reproducing and understanding LGM climate dating back more than 40 years (Alyea, 1972; Williams et al., 1973; Kutzbach and Guetter, 1986; Rind, 1987; Gates, 1976; Manabe and Hahn, 1977).

The complexity of processes and feedbacks that directly or indirectly influenced the LGM climate makes the task of reproducing it via dynamical models particularly challenging. Over the years, this has led to models of increasing complexity being applied to the study of the LGM climate, from intermediate complexity models to fully coupled Earth system models (ESMs) (Braconnot et al., 2012; Harrison et al., 2014, 2015; Annan and Hargreaves, 2015; Masson-Delmotte et al., 2006; Izumi et al., 2013; Li and Morrill, 2013; Lambert et al., 2013; Schmidt et al., 2014; Otto-Bliesner et al., 2007; Muglia and Schmittner, 2021, 2015; Menviel et al., 2017; Sime et al., 2013; Kageyama et al., 2013; Hargreaves et al., 2011). However, these increases in model complexity have not generally led to improved model performance when compared against proxy data (Harrison et al., 2014; Annan and Hargreaves, 2015). In particular, the latest results of the Paleoclimate Modelling Intercomparison Project phase 4 (PMIP4) for the LGM show that regional climate biases that characterized previous generation of climate models, such as underestimation of winter extratropical cooling and precipitation changes, are still present in the more recent model developments (Kageyama et al., 2021).

In recent years, regional climate models (RCMs) have been employed for the study of the LGM climate, mainly motivated by their improved representation of local processes and having a spatial resolution better matched to those of proxy reconstructions with respect to general circulation models (GCMs). In particular, several studies have shown that RCMs for the LGM can improve the simulated climate in comparison with the driving GCM simulation (Strandberg et al., 2011; Ludwig et al., 2019). Recently, a series of LGM studies have revealed that areas with mountainous terrain, such as the Alps, profit from the use of an RCM with convection-permitting resolution (Velasquez et al., 2020, 2021, 2022). They have also highlighted the important role of land-surface characterization for the representation of the LGM climate over Europe.

It is important to acknowledge that the added value of RCMs for paleo- applications is still under debate. While some studies highlight that an RCM signal is too dictated by its driving GCM, playing a major role (Armstrong et al., 2019), others consider the relatively low computational demands of RCMs as an added value compared to their driving GCMs, since it allows for a more comprehensive characterization of uncertainties related to changes in soil and surface features (Russo et al., 2022). Both points are actually very important for the LGM, for which high uncertainties in RCM simulations may derive from the “garbage in–garbage out” effect related to the imposed GCM signal and from differences in the characterization of surface features such as land cover and ice height (Kjellström et al., 2010; Strandberg et al., 2011). These sources of uncertainty have to be properly taken into account when willing to assess the climate of the LGM using an RCM.

In this study, we present a series of LGM sensitivity experiments performed over Europe with the RCM WRF 3.8.1 (Skamarock et al., 2008; Powers et al., 2017) but including some important technical developments with respect to the same model version used in the former studies of Velasquez et al. (2020, 2021, 2022). These developments are crucial for ensuring the model's accurate application in paleoclimate studies. The performed experiments are built considering two different land-cover datasets, as well as changes in continental and Alpine glaciers extent for both WRF and its driving GCM. Model results are evaluated against the newly developed pollen-based reconstruction database for the European LGM climate of Davis et al. (2022). The main goals of the study are the evaluation of a new model version against proxy reconstructions, the characterization of model uncertainties resulting from changes in the simulations setup relative to land cover and ice height and the assessment of the role of convection-permitting simulations for paleoclimate studies. At the same time, with a main focus on the Alpine LGM climate, the study is aimed at producing high-resolution outputs to force a glacier model and reconstructing the extent of glaciers during glacial maxima (Jouvet et al., 2023).

Section 2 gives a general overview of the different methods and data used in this study. In Sect. 3, the results are presented: first a comparison against the proxy-based reconstructions is conducted for summer and winter values of temperatures and precipitation in Sect. 3.1; secondly, the role of different uncertainties is considered for the same seasons and variables in Sect. 3.2; thirdly, some considerations on the role of convection-permitting resolutions are drawn in Sect. 3.3. Finally, a summary of the obtained results is discussed in Sect. 4.

2 Model, data, and methods

The results of this study are based on a set of simulations performed with the RCM WRF 3.8.1 (Skamarock et al., 2008; Powers et al., 2017). In the following, first the driving data used to conduct the performed simulations are introduced in Sect. 2.1. Successively, the general model setup is presented in Sect. 2.2, highlighting the main differences and improvements applied for the use of WRF in paleoclimate studies with respect to the same model version employed by Velasquez et al. (2020, 2021, 2022). Then, the complete set of performed LGM experiments is described in detail in Sect. 2.3. Finally, the new proxy reconstructions used for the assessment of the LGM model performance are introduced in Sect. 2.4.

2.1 Driving data

RCMs need climate information at their lateral boundaries. Here, outputs of a series of equilibrium climate simulations performed in a model chain of the Community Earth System Model (CESM 1.2 version, Hurrell et al., 2013) are used. CESM is a state-of-the-art Earth system model (ESM) consisting of different model components for the description of the atmosphere, land surface, ocean, and sea ice. In a first step, a fully coupled simulation at a horizontal resolution of 1.9° × 2.5° for the atmospheric and the land components and of nominal 1° × 1° for the ocean and the sea-ice components are executed (Buzan et al., 2023). Then, time-varying sea-ice mask and sea surface temperatures are derived from the simulation and are prescribed, in a successive step, to another experiment using only the atmosphere and land components of the CESM over an integration time of 24 years, at a horizontal resolution of 1.25° × 0.9°. Through these two steps, it is possible to reach a high spatial resolution, allowing for a realistic representation of the large-scale atmospheric circulation (Merz et al., 2015) while keeping the computational costs low.

For all of the CESM LGM simulations employed in this study, values of the orbital parameters (obliquity 22.949°, eccentricity 0.018994, and precession 294.425°), greenhouse gas concentrations (CO₂ 190 ppm, CH₄ 375 ppb, N₂O 200 ppb), and land use changes for the corresponding period at 21 ka are used following the directives of the Paleoclimate Model Intercomparison Project 4 (PMIP4, Kageyama et al., 2017).

The reference simulation is performed using the ICE-6G ice-sheet reconstruction of Peltier et al. (2015), the modification of which follows the setup for the LGM PIMP4 protocol (Zhu et al., 2017). Moreover, three additional sensitivity experiments are performed (Buzan et al., 2023), where only the ice-sheet height of the main Northern Hemisphere ice sheets (Laurentide, Greenland, Fennoscandia) is varied, respectively, to 33 %, 67 %, and 125 % of their height as derived from the dataset of Peltier et al. (2015). Additional details on the CESM model set up are presented in Buzan et al. (2023).

Figure 1 Maps of the topography and the different nested simulation domains.

2.2 WRF general model setup

The setup of the performed WRF simulations is partly based on the same model setup of earlier studies using the same model version 3.8.1 (Velasquez et al., 2020, 2021, 2022). This setup considers four two-way nested domains, with horizontal resolutions of 54, 18, 6, and 2 km, respectively (Fig. 1). The domains cover the entire European region, with a focus on the Alps. For the time integration, an adaptive time step is used with a minimum time step of 108, 36, 12, and 4 s, respectively, for the domains going from coarser to higher resolutions. The simulations are conducted considering a total of 40 vertical eta levels in the atmosphere, and four layers in the soil, with varying vertical resolution. The Kain–Fritsch cumulus convection scheme (Kain and Fritsch, 1990, 1993; Kain, 2004) is used for the coarser domains d01 and d02 of Fig. 1. In the inner domains d03 and d04, the convection parameterization is switched off since higher resolutions permit obtaining an explicit representation of convective processes. For all four domains, the land surface model NOAH-MP (Yang et al., 2011a) is used. For long-wave and short-wave radiation, the Dudhia and the RRTM schemes are respectively employed. The different LGM GHG values cannot be set in the model namelist but are passed directly to the corresponding RRTM shortwave radiation module prior to the model compilation, consistently with the values of the CESM simulations, as specified in Table 1. A general overview of the model setup employed for the different domains is presented in Table 2.

Table 1 Orbital parameters and greenhouse gas concentrations in PI and LGM periods.

Sakaguchi and Zeng (2009)

Table 2 General description of the WRF model setup of conducted LGM simulations.

Having introduced the main features of the WRF simulations in common with Velasquez et al. (2021), we successively describe the main differences with respect to this former study. The initial setup of Velasquez et al. (2021) included a glacier scheme of NOAH-MP, which allowed ice-phase changes in the soil, that was found to be incorrect, producing unrealistic soil temperatures. The here presented simulations use a slab ice scheme of NOAH-MP that does not consider ice-phase changes in the soil, resulting in more realistic soil temperatures than the former option (not shown). Another major improvement with respect to the model version of Velasquez et al. (2021) considered in this study, particularly relevant for the application of WRF to paleoclimate studies, concerns the model representation of changes in the orbital parameters of the Earth on millennial timescales. In fact, while different values of the obliquity of the orbit can be passed directly into the code of the default model version 3.8.1, this is not possible for the different values of the parameters representing changes in the eccentricity of the Earth's orbit around the Sun and in the precession of the equinoxes. For this reason, a FORTRAN subroutine already employed in WRF version 4.1.2 in Ludwig and Hochman (2022) is implemented in the main radiation module of the model to include the full orbital forcing of the LGM experiments. The routine allows for scaling the value of the solar constant depending on the effective position of the Earth on its orbit, taking into account changes in the orbital parameters. It has been developed from a subroutine used in several other paleoclimate studies with the RCM COSMO-CLM (Fallah et al., 2016, 2018; Prömmel et al., 2013; Russo and Cubasch, 2016; Russo et al., 2022) and is based on an original subroutine used in the GCM ECHAM5 (Roeckner et al., 2003), considering basic Kepler laws only and using the long-term series expansions of Laskar et al. (1993). The effect of the implemented subroutine on the seasonal pattern of insolation is here assessed for the coarsest study domain d01. Figure S1 of the Supplement shows the differences in incoming radiation on top of the atmosphere calculated for each day of a year between the LGM and the PI periods, for the default (left) and the new (right) version of the model. The results with the new treatment of orbital parameters (Fig. S1, right) are now in line with the expected seasonal pattern of insolation at the LGM over the study domain (Kageyama et al., 2017), substantially different than the one obtained from the original model version (Fig. S1, left).

2.3 WRF sensitivity experiments

The different CESM simulations with a spatial resolution of 1.25° × 0.9° introduced in Sect. 2.1 are used as initial and boundary conditions to run WRF. A total of five simulations with convection-permitting resolutions are performed with the WRF model for the LGM, each one covering approximately 11 years. For each of these simulations, given the extremely high computational demands of very high-resolution, each period is divided into three sub-periods of different length, between 52 and 40 months and considering the first 4 months as spinup time. In the end, the results of these sub-periods are joined together to retrieve a mean climatology for each of the performed experiments, based on the mean calculated over a total of 10 years. The approach is similar to the one of previous studies such as the one of Velasquez et al. (2020, 2021, 2022). It might not allow for exhaustively assessing the interannual variability in a model over a given period of study, due to the short time frames considered for each experiment. Nevertheless, since the interest of this study is mainly climatological mean values, this represents a good compromise between high demand in computational resources and the length of a signal to be used for building up a climatology. In support of this choice, the results of a 31-year long simulation of Velasquez et al. (2021) are also considered.

Table 3 Description of conducted LGM sensitivity experiments with WRF.

To perform experiments with WRF under LGM conditions, further modifications to the surface boundary conditions are necessary. Changes to land cover, soil composition, ice extent, ice height, and the land–sea mask are applied here to the original present-day WRF datasets. A reference simulation, here referred to as DEF (Table 3), is performed using the outputs from the reference CESM simulation of Buzan et al. (2023) as boundary conditions. For the DEF experiment, the LGM topography is prescribed directly from the dataset of Peltier et al. (2015). In particular, the maximum ice height in the period from 24 and 18 ka as derived from Peltier et al. (2015) is considered in this case. For the Alpine region, in the DEF experiment the ice-cap elevation is derived using the Parallel Ice Sheet Model (PISM, Khroulev and the PISM Authors, 2020) forced by the WRF LGM simulation of Velasquez et al. (2021), and tuned such that the glacial maximum area fits the reconstructed one by Ehlers et al. (2011). The land–sea mask of the DEF simulation is derived from Peltier et al. (2015), also considering the largest extension of ice in between 24 and 18 ka.

Figure 2 Orographic differences in the simulations (from top to bottom) ICE125, ICE67, and ICE33, calculated for the coarser domain of study d01 with respect to the reference simulation DEF. The different contours are drawn at intervals of 300 m. Dotted red contours represent negative values, while continuous green contours represent positive values.

Three additional simulations are conducted using in each case the different boundary data derived from the CESM simulations with perturbed LGM ice height over the higher latitudes of the Northern Hemisphere (see Sect. 2.1). Consistently with the driving data, for these three simulations the height of the northern latitudes ice sheets of the DEF simulation over the considered domain of study is increased by 25 %, and reduced by 67 % and 33 %. These simulations will be referred to, in the following text, respectively, as 125ICE, 67ICE, and 33ICE (Table 3). The orographic differences in each of these simulations with modified ice sheet, calculated for domain d01 with respect to the DEF experiment, are presented in Fig. 2. The most pronounced differences in prescribed topography are evident over the inner part of Scandinavia, with values lower than −1800 m obtained from the comparison between experiments 33ICE and DEF. In all the different cases, the ice-sheet extent (represented as white shadings in Fig. 3) is kept constant.

Figure 3 LGM land cover and biome classes for the coarser domain of study d01, passed as input to the DEF (a) and BIOME (b) simulations. Besides land ice cover indicated as ICE, all the legend acronyms of the different biome species derived from LPX-Bern are specified in Table 4.

Table 4 List of biomes from LPX-Bern (left) and the corresponding ones from MODIS they have been converted to (right).

For the LGM land cover, the LGM simulation by Velasquez et al. (2021) is used offline to drive the dynamic vegetation Land surface Processes and eXchanges model of the University of Bern (LPX-Bern v1.4) (Lienert and Joos, 2018; Spahni et al., 2013; Stocker et al., 2013). LPX-Bern features dynamic carbon, water, and nitrogen cycles and vegetation is internally represented by 10 plant functional types (PFTs; eight tree and two herbaceous) competing for resources and adhering to bioclimatic limits. Input variables for LPX-Bern are the fraction of land in a grid cell and monthly values of total precipitation, near-surface temperature, and net downward shortwave radiation at surface. Furthermore, for nitrogen deposition preindustrial values of the forcing product (Tian et al., 2018) closest to the model grid cell are used. LPX-Bern horizontal resolution is adapted to match the horizontal resolution of the respective WRF domains. For using the LGM maps of vegetation derived from LPX-Bern into WRF, the 12 biomes of LPX-Bern need to be translated into one of the 16 classes of WRF (based on the data from the US Geological Survey (USGS)). The selected BIOME correspondence between the two different datasets is presented in Table 3.

All the four simulations introduced above (namely DEF, ICE33, ICE67, and ICE125) use the same vegetation map derived with LPX-Bern from the simulation of Velasquez et al. (2021). Additionally, in order to acknowledge the possible effect of uncertainties in the prescribed land cover, a simulation is performed by driving offline LPX-Bern with the climatological outputs of the DEF simulation. This simulation will be referred to as the BIOME simulation. The DEF and the BIOME experiments are run using the same CESM boundary conditions at 1° × 1° resolution. A map of the two different land cover datasets used for the conducted experiments is presented in Fig. 3. For a large part of the domain there are noticeable differences in vegetation, both on a continental and local scale. The two datasets seem to correspond only over a remote part of northeastern Europe. For this reason, the designed experiments are likely to provide important information on the role of uncertainty related to the use of a different land cover in RCM simulations of European climate at the LGM. Present-day soil categories are used for all the LGM simulations, with some modifications: all the soil layers corresponding to a point covered by glaciers are set to ice. Consistently, for a given grid box, the same soil category of the upper soil layer is used for all soil levels at different depths.

In order to conduct a more quantitative evaluation of WRF against the pollen-based reconstructions taking into account different model uncertainties, the five experiments introduced above (i.e., DEF, ICE33, ICE67, ICE125, and BIOME) are considered together as an ensemble of simulations. Here, as model uncertainty, we consider the range of outputs that is obtained using the same model but applying changes to the model chain setup inherent to land cover and ice height.

Finally, in addition to the already described simulations, an experiment with the same setup of the DEF simulation but using only two nested domains, D01 and D02, down to a spatial resolution of 18 km and with the convection parameterization switched on (i.e., convection is not explicitly solved), is performed. This experiment is indicated as DEF_parconv in Table 3 and is conducted with the goal of better assessing the role of convection-permitting simulations for the representation of the European LGM climate.

2.4 Proxy reconstructions dataset

The model simulations are evaluated against a new pollen-based reconstruction of LGM seasonal and annual temperature and precipitation by Davis et al. (2022). This reconstruction is based on pollen data from 63 pollen records that were selected for the quality of their dating control over the LGM period (21 ka cal ±2 kyr). Records were excluded if they fell within category 7 (“poor”) of the 1–7 COHMAP scale of chronological quality based on the proximity of dates to the LGM time slice. Only absolute dating methods (e.g., radiometric dating such as radiocarbon) were considered, and relative dates based on pollen correlation were excluded. The record selection process, therefore, excluded many of the pollen records used in previous studies, reflecting the generally poor quality of the dating control in LGM pollen records, including 10 of the 18 records used in the PMIP benchmarking dataset (Bartlein et al., 2011). The reconstruction itself uses a standard modern analogue technique (MAT) that finds the best six analogues amongst a dataset of over 8000 modern pollen samples from the Eurasian Modern Pollen Database (Davis et al., 2020). The MAT method has been used for previous LGM pollen–climate reconstructions (Peyron et al., 1998) but it has not been known how much the method may be influenced during the LGM by problems such as having no modern analogue vegetation (and climate), as well as changing plant physiological responses to climate associated with lower atmospheric CO₂. For this reason, the method was evaluated against a previous pollen–climate reconstruction based on inverse modeling (Wu et al., 2007), which is designed to account for many of these problems, as well as a chironomid-based reconstruction (Samartin et al., 2016). The result of this evaluation showed little difference between the MAT reconstruction and these other methods, indicating that the MAT method is reliable for this time period and location. This can be attributed in part to advances in the size and spatial coverage of the modern pollen training set, which is an order of magnitude larger than those used in previous MAT reconstructions, permitting good vegetation analogues (chord distance <0.3, Huntley, 1990) to be found for all fossil pollen samples. Importantly, the uncertainty estimates of the considered pollen-based reconstructions are considered in their comparison against the results of the presented model simulations.

Figure 4 LGM winter (a, c) and summer (b, d) climatological values of mean near-surface temperature (a, b) and precipitation (c, d) as derived from the reference WRF LGM simulation, for the domain d01 of Fig. 1, at coarser resolution. The circles represent the values as derived from pollen-based reconstructions.

3 Results and discussion

3.1 Comparison of the reference simulation against pollen-based reconstructions

In Fig. 4 we present the mean climatologies of winter and summer temperatures and precipitation of the DEF LGM simulation over the coarsest domain of the study, together with the corresponding values derived from the pollen-based reconstructions. Winter temperatures below 0 °C characterize almost the entire domain for the DEF experiment, except for a few points over southern Spain and northwestern Africa. Particularly cold simulated temperatures are evident in winter over the northeastern part of Europe, with values generally below −25 °C. Temperatures remain lower than 0 °C in the model in summer over north and northeastern Europe, over areas largely covered by glaciers. Simulated temperatures rarely exceed a mean value of 20 °C in summer, only over specific areas of southern Europe and northern Africa. The values of the pollen-based reconstructions, displayed in Fig. 4 as filled circles, show a similar behavior, even though their coverage is mainly limited to central and southern Europe. In winter, the pollen data are also characterized by a northeast-to-southwest gradient in temperatures, with values consistent with the ones of the WRF DEF experiment. In summer, reconstructed temperatures are similar to the ones simulated by WRF for the entire domain, presenting in particular a similar north-to-south gradient. The simulated precipitation pattern in winter is characterized by high values in the western areas of the domain (Fig. 4, lower row), for which moisture advection from the Atlantic Ocean plays a major role (Beghin et al., 2016; Ulbrich et al., 1999), up to seasonal values of 500 mm/season in some cases. In summer, drier conditions are evident over almost the entire domain, except for few areas with complex topography such as part of the Alps and western Scandinavia, where precipitation exceeds 300 mm/season. The pollen-based reconstructions show a slightly different picture compared to the DEF simulation for precipitation in winter. While the western coastal areas are wetter than eastern Europe according to the pollen, pronounced dry conditions are evident over the central part of the domain, particularly at the northern edges of the Alps. Additionally, the range of values of winter precipitation derived from the pollen-based reconstructions is generally narrower than in the case of the DEF simulation, with maximum values rarely exceeding 250 mm/season. For summer precipitation, on the other hand, there is a better agreement between the simulated data and proxy reconstructions than in winter, except for specific areas, such as again around the Alps. In summer, the simulated values of precipitation are in the same range as the climate reconstructions.

Figure 5 Maps for the bias of winter (a, c) and summer (b, d) temperatures calculated between the WRF reference simulation and the pollen-based reconstructions. Upper row (a, b) shows the results for the coarser-resolution domain d01, while the bottom row (c, d) presents the same analysis for the higher-resolution domain d04. The colored circles show the biases calculated between the two datasets. In the case of the bias being smaller than 1 standard deviation of the reconstructions, a black dot is added at the center of the circle. Considering then the five-member ensemble, in the case of the bias of at least one of these members against the reconstructions being smaller than 1 standard deviation of the latter, then the circle is highlighted in red. Finally, the gray scale represents the maximum differences for each point of the domains, between the different members of the ensemble.

Figure 6 Maps for the bias of winter (a, c) and summer (b, d) seasonal precipitation calculated between the WRF reference simulation and the pollen-based reconstructions. Upper row (a, b) shows the results for the coarser-resolution domain d01, while the bottom row (c, d) presents the same analysis for the higher-resolution domain d04. The colored circles show the biases calculated between the two datasets. In the case of the bias being smaller than 1 standard deviation of the reconstructions, a black dot is added at the center of the circle. Considering then the five-member ensemble, in the case of the bias of at least one of these members against the reconstructions being smaller than 1 standard deviation of the latter, then the circle is highlighted in red. Finally, the gray scale represents the maximum differences for each point of the domains, between the different members of the ensemble.

3.2 Consideration of different model uncertainties

A more quantitative assessment of the differences in seasonal values of temperature and precipitation between the DEF WRF simulation and the pollen-based reconstructions is provided in Figs. 5 and 6, respectively. In these figures, the temperature and precipitation derived from the additional WRF simulations BIOME, 33ICE, 67ICE, and 125ICE of Table 3 (forming, together with the DEF simulation, a five-member ensemble) are also considered. The comparison between the pollen-based reconstructions and the model outputs is additionally conducted on both the coarse-resolution domain, and the inner domain d04 of Fig. 1 in this case. A circle with the values of the differences between the reference simulation DEF and the pollen-based reconstructions is drawn for the proxy locations. If the bias of the reference simulation DEF lies within 1 standard deviation of the uncertainty values of the pollen-based reconstructions, then a black dot is added at the center of a given circle. Successively, if at least one of the five-member ensemble lies within 1 standard deviation of the pollen values, then the circle is highlighted in red. Finally, the maximum absolute differences in the considered variables calculated between the different ensemble members are plotted for each point of the domain on a gray scale in order to provide an assessment of the ensemble spread and model uncertainty.

Figure 5 confirms that, for both winter and summer, there is a very good agreement in temperatures between the performed simulations and the values derived from the proxies. The consideration of the different simulations with the applied changes in the model configuration allow for obtaining results that are very close to the evidence of the proxies for almost the entire points of the domain, in both seasons and for both considered variables. For winter temperature, the bias of the reference (DEF) exceeds in very few cases an absolute value of 5 °C. Interestingly, despite the large changes in the configuration of the different ensemble members, their results in terms of winter temperature are relatively close for almost the entire domain of study (see also Fig. S2). Differences between the different ensemble members larger than 4 °C are evident for winter temperature mainly over the northeastern part of the domain and Iceland, with values of the different ensemble members spreading here over a range of more than 10 °C. For summer, again, the mean model temperatures are very close to the values of the pollen-based reconstructions over the entire domain (Fig. 5, right column). The DEF simulation, in particular, shows a small bias against the proxies over the entire domain, with values rarely exceeding 4 °C and well within the range of the pollen uncertainty. The main exceptions are some points over the area of present-day France, for which a closer agreement is reached only when considering other members of the ensemble. In summer, more notable and spatially heterogeneous differences in temperatures are evident between the different ensemble members than in winter. In particular, large differences, of up to 14 °C, are found between the different experiments for areas at the border of glaciers over Scandinavia and Iceland, mainly due to a reduction in the height of glaciers over these regions (Fig. S3).

Analyses for the inner Alpine domain (Fig. 5, bottom row) confirm an overall good match between temperatures derived from the performed simulations and the pollen-based reconstructions, for both summer and winter. In this case, the main differences between the different ensemble members are evident in winter over the Po valley, exceeding in some cases 6 °C. Just two points of the pollen reconstructions, in summer, have temperature values outside the range of the ensemble members for the inner domain of study. This gives high confidence on the use of the presented results for the study of the LGM climate of the area.

For precipitation (Fig. 6), the agreement between model simulations and proxy reconstructions is not as good as in the case of temperature. In winter, the model seems to be generally too wet, in particular over the northern Iberian Peninsula, present-day France, the Alps, and continental areas of present-day Greece and Turkey, with seasonal values of precipitation exceeding the ones derived from pollen-based reconstructions by more than 300 mm/season. For these areas, none of the performed simulations gets significantly closer to the range of values of the proxies, also when taking into account pollen uncertainties. It is important to note here that most of the locations for which the model presents the largest positive biases in winter precipitation are at the boarder of mountainous areas covered with glaciers, such as the Alps and part of the Pyrenees. This positive bias could partly be related to uncertainties in the extent or height of the prescribed glaciers, not properly resolved at the scales characteristic of the considered pollen data. This is confirmed by the fact that the bias for the two points northwest of the Alps is reduced when considering the model results over the Alpine region at an horizontal resolution of 2 km (Fig. 6, bottom row). In winter, the areas for which the DEF simulation presents higher values of precipitation, i.e., western Europe (Fig. 4, lower row), are the same as those for which higher differences in the ensemble members are evident, exceeding in most of these areas 150 mm/season. These pronounced differences between the different ensemble members are mainly dictated by the experiments with changes in the height of glaciers (Fig. S4). A similar overlapping of the largest differences in the ensemble members over regions presenting the highest values of precipitation in the DEF simulation is also evident in summer. In this case though, there is a better agreement between the results of the performed simulations and the pollen-based reconstructions, although the model seems to be generally drier, in particular over the Balkans, where the differences against the pollen-based reconstructions are in some cases below −200 mm/season. Notable in this case is that changes in land cover seem to be very important for the representation of summer precipitation over the Mediterranean region, in particular over the Balkans (Fig. S5), where the BIOME simulation produces significant changes with respect to the DEF experiment.

For the inner domain of study there is also a relatively good agreement between the proposed simulations and the pollen-based reconstructions in terms of precipitation, for both seasons, but this is again worse than in the case of temperature (Fig. 5, lower row), at least when considering the number of points for which the simulations lie within the pollen uncertainty. Three important points need to be mentioned in the case of precipitation for the inner domain of study. Firstly, the differences between the different ensemble members are considerable over all of the Alpine range, in both winter and summer, with values extensively exceeding 200 mm/season. This confirms again the relevance of the considered sources of uncertainty in the model for the investigation of the LGM climate of the region. Secondly, even though the simulations are out of the range of the considered pollen-based reconstructions over some points in one season, they get closer in the other. This indicates that some of the results, even if not perfectly matching the values of the proxies, are not completely off. Thirdly, a point also valid for the rest of the domain is that the pollen reconstructions are trained with present-day observations, and studies show that even observations nowadays tend to underestimate precipitation in particular areas characterized by complex terrain (Fu, 1996; Frei and Schär, 1998; Frei et al., 2003).

Finally, to confirm that the differences between the different ensemble members are robust and not simply an effect of the selection of a relatively short simulation period of 10 years, from which the climatological values are obtained, we have also calculated the maximum differences in the climatological values of 20 different periods of 10 years derived from the 31-year long simulation of Velasquez et al. (2021), for summer temperatures. The results, as illustrated in Fig. S6, indicate that the spread of this set of experiments is negligible compared to the one obtained when comparing the different realizations of our five-member ensemble.

Figure 7 Bias in winter (a, c) and summer (b, d) values of temperatures (a, b) and precipitation (c, d) calculated for the coarser-resolution domain d01, between the reference simulation and the one with only two nested domains, down to 18 km. The dots in each of the panels represent the points for which the values differ significantly, at a significance level of 0.05, according to the results of a KS significance test.

3.3 Considerations on the role of convection-resolving resolutions

In Fig. 7 the differences in the mean seasonal values of temperature and precipitation between the reference DEF simulation and the DEF_parconv experiment are calculated for each grid box of the outer domain d01. The main goal of such an analysis is to assess whether the explicit representation of convection over the inner part of the domain has a relevant impact on the simulated variables at both the local and continental scale, with respect to a simulation exclusively using a parameterized representation of convection over all the domain. In Fig. 7, for each grid-box where the daily mean anomalies of the considered variable and season can be considered significantly different according to the results of a Kolmogorov–Smirnov (KS) significance test at a significance level of 0.05, a dot is plotted.

Results show that for the Alpine region, the effects of higher resolution and the explicit representation of convection are significant, for both variables and seasons. For precipitation, in particular, the effects of the convection-permitting simulation are seen not only over the inner domains but also over the entire domain of study, presenting significant differences between the two simulations. These differences are in some cases of the same order of magnitude as the differences between the different sensitivity tests conducted by perturbing the glaciers ice height and extension (see Supplement, Fig. S2) and considering considerably different boundaries discussed in Sect. 3.2. Therefore, the consideration of convection-permitting resolutions for the simulation of European LGM climate has a significant impact on the results, comparable to the effect of other relevant sources of uncertainty, not only at a local scale but over the entire European domain.

4 Conclusions

In this study, we conduct an evaluation of a version of the RCM WRF 3.8.1, including a series of important technical developments necessary for paleoclimate applications, against newly available pollen-based reconstructions for the LGM climate of Europe and the Alpine region. Through the performance of a set of sensitivity experiments taking into account the role of different large-scale and surface model error sources, such as changes in the imposed northern hemispheric ice height and land cover, we aim to assess the general performance of the model. At the same time, we quantify the possible effect of changes in the model setup on the obtained results, highlighting where outcomes of RCMs can be considered more robust and where factors such as error in the representation of surface features could play a major role in the reconstruction of the European LGM climate.

Results show that the model ensemble produced by changing surface drivers and using different boundaries is in very good agreement with the considered pollen-based reconstructions for summer and winter values of temperature. For precipitation, the agreement is also relatively good, but in this case none of the produced ensemble members gets as close to the proxies over the entire domain of study as in the case of temperature. In particular, the model results are too wet in winter over western Europe and too dry in summer over the eastern part of Europe, with respect to the pollen-based reconstructions, even when considering the uncertainty associated with the considered pollen-based reconstructions. The model shows a positive northeast-to-southwest gradient of temperatures in winter over Europe at the LGM, in accordance with the picture drawn from the pollen-based reconstructions. Most of the European domain presents mean winter values below 0 °C as derived from the model simulations. A positive north-to-south gradient characterizes simulated summer temperatures, again similarly to evidence from the proxies (Russo, 2023). In summer, mean simulated temperatures rarely exceed 20 °C. The LGM WRF precipitation shows particularly high values, exceeding in some cases 400 mm/season over Atlantic regions of western Europe and the Alps. In summer, the pattern of precipitation in the model is quite homogeneously distributed, with values rarely exceeding 150 mm/season, except for specific areas characterized by complex topography, such as Scandinavia and the Alps. The consideration of different northern hemispheric ice height and land cover types in the model leads to largely differing precipitation values over specific areas, with differences between the different experiments exceeding 200 mm/season in some cases. This allows for increasing the match between the model and the proxies, at the same time highlighting the importance of the considered sources of uncertainty for the simulation of European climate at the LGM over different regions. Errors associated with changes in the ice height of the higher latitudes of the Northern Hemisphere produce very pronounced changes in both winter and summer precipitation over western Europe and the Alpine region. Changes in land cover also play an important role in summer precipitation over the Mediterranean basin. The considered changes have a relatively limited impact on winter temperatures over almost the entire domain, playing instead a larger role in summer temperatures, in particular over areas at the boarders of glaciers, with differences between the different ensemble members reaching 14 °C in some cases. Finally, the use of four two-way nested domains, enabling reaching a spatial resolution in the inner domain of approximately 2×2 km and explicitly representing convection, leads to important changes in model results with respect to an experiment with parameterized convection over the inner part of the domain. Significant differences between the two experiments are evident not only over the inner Alpine region but also over the outer domain of study. These changes are particularly pronounced for precipitation rather than for temperature, being in the range of the changes relative to the use of different ice height and land cover in this case.

The paper sets the basis for an improved representation of the LGM climate of Europe. It introduces some important technical developments useful for the application of WRF to the study of the LGM with respect to previously employed model versions, allowing for obtaining seasonal values of temperature and precipitation in very good agreement with newly available pollen-based reconstructions. Additionally, it permits assessment of the role of different sources of model uncertainty, as well as where and for which variable largest differences in model results might be expected as a result of a different model setup. Finally, the study enables assessment of the effect of convection-permitting resolutions for the study of the LGM climate at both the local and continental scale.