Applying an isotope-enabled regional climate model over the Greenland ice sheet: effect of spatial resolution on model bias

In order to investigate the impact of spatial resolution on the discrepancy between simulated δ¹⁸O and observed δ¹⁸O in Greenland ice cores, regional climate simulations are performed with the isotope-enabled regional climate model (RCM) COSMO_iso. For this purpose, isotope-enabled general circulation model (GCM) simulations with the ECHAM5-wiso general circulation model (GCM) under present-day conditions and the MPI-ESM-wiso GCM under mid-Holocene conditions are dynamically downscaled with COSMO_iso for the Arctic region. The capability of COSMO_iso to reproduce observed isotopic ratios in Greenland ice cores for these two periods is investigated by comparing the simulation results to measured δ¹⁸O ratios from snow pit samples, Global Network of Isotopes in Precipitation (GNIP) stations and ice cores. To our knowledge, this is the first time that a mid-Holocene isotope-enabled RCM simulation is performed for the Arctic region.

Under present-day conditions, a dynamical downscaling of ECHAM5-wiso ($\mathrm{1.1}{}^{\circ }\times \mathrm{1.1}{}^{\circ }$) with COSMO_iso to a spatial resolution of 50 km improves the agreement with the measured δ¹⁸O ratios for 14 of 19 observational data sets. A further increase in the spatial resolution to 7 km does not yield substantial improvements except for the coastal areas with its complex terrain. For the mid-Holocene, a fully coupled MPI-ESM-wiso time slice simulation is downscaled with COSMO_iso to a spatial resolution of 50 km. In the mid-Holocene, MPI-ESM-wiso already agrees well with observations in Greenland and a downscaling with COSMO_iso does not further improve the model–data agreement. Despite this lack of improvement in model biases, the study shows that in both periods, observed δ¹⁸O values at measurement sites constitute isotope ratios which are mainly within the subgrid-scale variability of the global ECHAM5-wiso and MPI-ESM-wiso simulation results. The correct δ¹⁸O ratios are consequently not resolved in the GCM simulation results and need to be extracted by a refinement with an RCM. In this context, the RCM simulations provide a spatial δ¹⁸O distribution by which the effects of local uncertainties can be taken into account in the comparison between point measurements and model outputs. Thus, an isotope-enabled GCM–RCM model chain with realistically implemented fractionating processes constitutes a useful supplement to reconstruct regional paleo-climate conditions during the mid-Holocene in Greenland. Such model chains might also be applied to reveal the full potential of GCMs in other regions and climate periods, in which large deviations relative to observed isotope ratios are simulated.

1 Introduction

Stable isotopes of water (HD¹⁶O and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) are fractionated during any phase transition. This fractionating process depends on temperature (Dansgaard, 1953; Craig and Gordon, 1965; Jouzel and Merlivat, 1984), and thus water isotopic ratios (expressed here in the usual δ notation, δD and δ¹⁸O with respect to the Vienna Standard Mean Ocean Water, V-SMOW) reflect the atmospheric conditions under which the fractionating process took place (Dansgaard, 1964; Merlivat and Jouzel, 1979; Gat, 1996). This process is generally utilized to reconstruct paleo-climate conditions such as past temperature changes, using isotopic ratios stored in climate archives (Dansgaard et al., 1969; Masson-Delmotte et al., 2005; Jouzel, 2013).

In Arctic regions like Greenland, ice cores constitute an exceptional climate archive. Over thousands of years, accumulated snow was solidified to ice, preserving the water isotopic ratios since the last interglacial period at some locations. Climate reconstructions based on these ice cores show that the climate conditions changed considerably in Greenland during the Holocene (here defined as the period between present day and 12 ka; Marcott et al., 2013). Between the early Holocene and the Holocene Thermal Maximum in the mid-Holocene (6 ka), a pronounced warm phase took place. Since then, temperatures steadily decreased until the late Holocene (Marcott et al., 2013; Moossen et al., 2015). In this context, the mid-Holocene is a period of particular interest, as by that time an Arctic warming had taken place due to orbital forcing variations and their related feedbacks on large-scale climate variations, which exhibits similarities to the strong recent Arctic warming. For Greenland, the mid-Holocene provides the opportunity to investigate the processes leading to this warming in more detail and to potentially obtain new insights about the future development of the Arctic region (Yoshimori and Suzuki, 2019).

While general circulation models (GCMs) are generally able to reproduce the direction and large-scale patterns of past climate changes (e.g., Timm and Timmermann, 2007; Smith and Gregory, 2012), they often fail to reproduce the magnitude of regional changes (Braconnot et al., 2012; Harrison et al., 2014) documented in various local climate archives. Thus, a scale gap might exist between the measured point information and the large-scale climate information generated by GCMs. The comparison of observational and GCM data can therefore be subject to considerable uncertainties (Felzer and Thompson, 2001).

Especially for structured landscapes, the spatial resolution in GCMs is often too coarse to resolve relevant local factors (Jost et al., 2005; Fischer and Jungclaus, 2011). Important properties like topography and surface conditions are consequently only represented in a generalized and imprecise form in climate simulations. In most cases, this does not adequately represent the complex characteristics of the land surface and its associated interactions with the atmosphere. For stable water isotopes, key physical processes of isotope fractionation are therefore not well resolved in coarse resolution GCMs, leading to differences between simulated and observational isotope data, especially in complex terrains (Sturm et al., 2005; Werner et al., 2011). Isotope-enabled GCMs are consequently not able to quantitatively reproduce regional changes in isotope ratios (e.g., Risi et al., 2010), and the simulated isotope ratios with GCMs exhibit in many cases larger deviations relative to observed ratios than the results of corresponding regional climate model (RCM) simulations. For instance, Sturm et al. (2007) were able to reduce the bias of simulated isotope ratios in precipitation through a regional downscaling of an isotope-enabled GCM run in South America. Comparable results were achieved by Sjolte et al. (2011) for isotope-enabled RCM simulations in Greenland.

Therefore, in the present study, isotope-enabled GCM simulation results for the Arctic region are dynamically downscaled with an isotope-enabled RCM to a higher temporal and spatial resolution. By means of such regional simulations, the spatial and temporal variability of the isotopic ratios in the Arctic is potentially increased, accounting for the heterogeneity of local conditions at the different ice core locations and the associated uncertainties. In this way, the impact of highly resolved local conditions on the spatial and temporal variability of isotopic ratios is investigated, and the impact of such small-scale variability on the discrepancy between simulated and observed paleo-climate conditions in the Arctic region is examined.

To explore this, the isotope-enabled version of the RCM COSMO-CLM (Rockel et al., 2008), COSMO_iso (Pfahl et al., 2012; Christner et al., 2018), is used. In a first step, the general suitability of COSMO_iso to be used for isotope applications in Greenland is assessed. For this purpose, near-surface temperatures and precipitation amounts simulated with the standard COSMO version are compared with observations in the Arctic region. Subsequently, the capability of COSMO_iso to simulate realistic water isotopic ratios for Greenland is tested by downscaling a global present-day simulation with an isotope-enabled GCM for the Arctic region. The GCM and RCM results are then compared to measured water isotope ratios in precipitation and snow pit samples. Afterwards, the tested isotope-enabled COSMO_iso model system is used to downscale an isotope-enabled GCM simulation for a mid-Holocene time slice. The simulated isotopic ratios are evaluated against Greenland ice core data. Such a dynamical downscaling of global isotope simulations for Greenland under mid-Holocene conditions is performed for the first time in the framework of this study.

2 Methods

2.1 COSMO_iso

2.1.1 Model description

In this study, simulated stable water isotope concentrations of HD¹⁶O and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ with isotope-enabled GCMs (Sect. 2.1.2), are regionally downscaled with COSMO_iso (Pfahl et al., 2012), an isotope-enabled version of the numerical weather prediction model COSMO (Consortium for Small-scale Modeling; Baldauf et al., 2011) (version 4.18). For the purpose of long-term climate simulations, isotope routines of COSMO_iso were implemented in COSMO-CLM (Rockel et al., 2008), the climate version of COSMO. In this context, the δD and δ¹⁸O ratios in the soil water and the surface-layer snow are simulated with TERRA_iso version 1 (Dütsch, 2017; Christner et al., 2018), the isotope-enabled version of the multi-layer Land Surface Model TERRA-ML (Schrodin and Heise, 2001) in COSMO. In several studies, COSMO_iso and TERRA_iso were successfully employed for the simulation of isotopic ratios in the midlatitudes (Pfahl et al., 2012; Aemisegger et al., 2015; Christner et al., 2018). In the present study, the model system will be applied to the Arctic region. For this, some additional modifications regarding the treatment of snow and ice had to be implemented in the model.

Snow albedo

The surface albedo of fresh snow is increased from 0.7 to 0.8 to improve the model agreement with measured values of short-wave reflectance and 2 m temperature at stations from the Cooperative Institute for Research in Environmental Sciences at the University of Colorado Boulder (CIRES) in central Greenland (Karremann and Schädler, 2021).

Snow layer thickness

In the standard configuration of COSMO, the Greenland ice sheet is treated as a constant mass of ice that is covered by a single snow layer. But in this model structure, dynamical processes within the ice sheet (flow, basal melt) are not included. As a result, the depth of the snow layer is constantly increasing and thus so is its heat capacity. To avoid this spurious model behavior, the snow layer depth is limited to 5 cm in this study. Using this value, realistic diurnal cycles of the 2 m air temperature could be simulated.

Marine regions with sea ice cover

To be able to simulate reasonable fractionation processes for marine regions with sea ice cover, a snow layer is also implemented on top of the sea ice (e.g., as suggested in Bonne et al., 2019). The isotopic composition of this surface snow layer is in this case set to the isotopic composition of the most recent precipitation.

Fractionation at snow covered surfaces

Isotope fractionation during sublimation from a surface snow layer is poorly understood. Several different processes are suggested to be involved, which are not yet taken into account in state-of-the-art isotope-enabled models (see, e.g., discussion in Christner et al., 2017), such as non-fractionating layer-by-layer sublimation (e.g., Ambach et al., 1968), kinetic fractionation during sublimation into sub-saturated air, a diurnal cycle of sublimation combined with fractionating vapor deposition on the snow (e.g., Steen-Larsen et al., 2014) and fractionating melt water evaporation combined with recrystallization of residual melt water (Gurney and Lawrence, 2004). To approximate this complex interplay of different influencing factors, in this study an equilibrium fractionation during sublimation from surface layer snow and sea ice is assumed. However, the authors are aware that this is a simplified description of isotope fractionation during sublimation.

2.1.2 Model simulation setup

The capability of COSMO_iso to realistically reproduce the fractionating processes of stable water isotopes in Greenland is evaluated. For this, the nudged simulation outputs (standard and isotopic) from an isotope-enabled atmospheric model ECHAM5-wiso (Werner et al., 2011) simulation are dynamically downscaled with COSMO_iso for the whole Arctic region. The data from the same ECHAM5-wiso simulation have been already used as boundary conditions for COSMO_iso simulations over Europe by Christner et al. (2018). The simulation outputs from ECHAM5-wiso are at a T106 horizontal spatial resolution ($\mathrm{1.1}{}^{\circ }\times \mathrm{1.1}{}^{\circ }$) and on 31 atmospheric vertical levels. The dynamical fields were nudged every 6 h towards ERA-Interim reanalysis data (Dee et al., 2011). Monthly varying sea surface temperatures and sea ice cover were prescribed as lower boundaries over sea, also based on the ERA-Interim data. The simulation period is 2008–2014. In order to guarantee a physically consistent transition between the coarse model resolution of the GCM and the fine model resolution of the RCM at the boundaries of the regional model domain, the model resolution is increased stepwise. This procedure is called nesting. In a first nesting step, the spatial resolution of COSMO_iso is set to $\mathrm{0.44}{}^{\circ }\times \mathrm{0.44}{}^{\circ }$, corresponding to 50 km×50 km in rotated coordinates (COSMO_iso_50km). In a second nesting step, an additional COSMO_iso simulation with a spatial resolution of $\mathrm{0.0625}{}^{\circ }\times \mathrm{0.0625}{}^{\circ }$ (corresponding to about 7 km×7 km) for Greenland (COSMO_iso_7km) is nested in the COSMO_iso_50km simulation. This high-resolution simulation covers the year 2011. In the COSMO_iso runs, the horizontal wind fields above the 850 hPa level are spectrally nudged (von Storch et al., 2000) towards the reanalysis-based dynamical fields of ECHAM5-wiso. This method ensures that consistent atmospheric boundary conditions build the framework for the fractionating processes simulated in COSMO_iso. The model domains of both simulations is shown in Fig. 1.

Figure 1 COSMO_iso model domain of the 50 km (blue) and the 7 km (black) simulation.

The same model chain is applied to the mid-Holocene period. Atmospheric fields have been retrieved from a mid-Holocene simulation of the fully coupled isotope-enabled Max Planck Institute Earth System Model (MPI-ESM-wiso, Cauquoin et al., 2019), whose atmospheric component is ECHAM6-wiso. The major ECHAM6 model changes compared to ECHAM5 include an improved representation of radiative transfer in the solar part of the spectrum, an improved representation of surface albedo, a new aerosol climatology and an improved representation of the middle atmosphere (Stevens et al., 2013). The ocean component is the Max Planck Institute Ocean Model (MPIOM, Jungclaus et al., 2013). With COSMO_iso, only a representative time slice of 30 years is simulated for this climate period, as the regional COSMO_iso simulations are computationally very expensive. The greenhouse gas concentrations and the orbital parameters are adapted according to the Paleoclimate Modelling Intercomparison Project 4 experiment design (PMIP4, Kageyama et al., 2018). The model domain of the COSMO_iso simulations is identical to the present-day simulations.

2.2 Observations

The capability of the isotope-enabled regional climate model COSMO_iso to reproduce measured isotopic ratios in Greenland is evaluated by comparing the simulation results to observational data. The simulated isotopic composition in precipitation is assessed by comparing the model results in the Arctic region with observed monthly data from the Global Network of Isotopes in Precipitation (GNIP) of the International Atomic Energy Agency and the World Meteorology Organization (IAEA, 2018) over the period 2008–2014 (Table 1). Furthermore, simulated δ¹⁸O ratios are compared to snow pit samples collected during the North Greenland Traverse (Fischer et al., 1998; Weißbach et al., 2016a) and top core samples from five ice core locations (Renland, Vinther et al., 2008; Neem, Masson-Delmotte et al., 2015; GISP2, Grootes and Stuiver, 1997; Summit, Fischer, 2003; SE-Dome, Furukawa et al., 2017). The station numbers assigned to the respective samples within this study, as well as their locations and δ¹⁸O values are summarized in Table 2. Since all snow pit samples cover different time periods, the present-day δ¹⁸O values (black numbers in Table 2) are calculated as an average of all available δ¹⁸O values measured between 1940 and 2014. With this procedure, uncertainties in snow pit samples and top ice core samples, associated with post-depositional diffusion and wind erosion and the resulting constraints in analyzing annual and interannual top ice core data (e.g., Johnsen et al., 2000), can be neglected. However, further uncertainties in snow pit samples and ice core data remain regarding the timescale assignment (Steig et al., 2005) and the spatial variability (Weißbach et al., 2016b).

Table 1 List of GNIP stations used in this study.

Table 2 Description of the snow pit and ice core samples used in this study. The present-day δ¹⁸O values are calculated as an average of all available δ¹⁸O values measured in snow pit samples between 1940–2014. The mid-Holocene δ¹⁸O values are calculated as an average of the measured δ¹⁸O values in ice cores over the period 5.5–6.5 ka. Black numbers indicate present-day δ¹⁸O values and the respective standard deviations, and bolded numbers indicate mid-Holocene values and the corresponding standard deviations.

Since both snow pit samples and top core samples from ice cores represent an integrated signal of the isotopic composition in precipitation, we computed modeled annual mean δ¹⁸O values in precipitation and compared the multi-year 2008–2014 model mean to the observed values. For the calculation of the yearly modeled mean values, the modeled δ¹⁸O in precipitation is weighted with accumulation rate, i.e., months with high precipitation amounts get a higher weight compared to months with small precipitation amounts.

Ice core samples are also used to evaluate the simulated isotopic ratios for the mid-Holocene. Besides the already-mentioned Renland and GISP2 samples, two more ice core samples, namely GRIP and NGRIP (Vinther et al., 2006), are used for the model evaluation. The mid-Holocene δ¹⁸O values (blue numbers in Table 2) are calculated as an average of the measured δ¹⁸O values in ice cores over the period 5.5–6.5 ka.

3 Results

3.1 Present day

3.1.1 Standard climatological parameters

In a first step, the general capability of the COSMO model to reproduce observed standard climatological parameters in present-day simulations for Greenland is assessed. For this purpose, the results of an ERA-Interim reanalysis-driven (Dee et al., 2011) simulation with the standard COSMO model (without isotope application) are compared with observations collected by the Danish Meteorological Institute (DMI). Evaluated are the yearly mean 2 m temperatures (Fig. 2a) and the yearly mean precipitation sums (Fig. 2b) over the period 1995–2015. Both simulated 2 m temperatures and simulated precipitation amounts are in good agreement with the DMI observations. The simulated 2 m temperatures coincide especially well with the observed values. For the precipitation sums, the spread of simulated and observed values is higher than for the 2 m temperatures, a feature generally occurring in weather and climate simulations. Thus, the COSMO model is generally able to simulate reasonable climate conditions in Greenland and can therefore be used for isotope applications in this region. A detailed analysis of the standard COSMO model performance in the Arctic region is presented in Karremann and Schädler (2021).

Figure 2 Simulated yearly mean (a) 2 m temperatures and (b) precipitation sums of a standard COSMO simulation, driven with ERA-Interim, for Greenland over the period 1995–2015 compared to DMI observations. The solid black line is the 1:1 line in (a, b).

3.1.2 Comparison of simulated δ¹⁸O data to station data

Figure 3 shows the yearly mean δ¹⁸O values for the period 2008–2014 for Greenland, simulated with COSMO_iso_50km (a) and ECHAM5-wiso (b). Additionally, the locations and the observed δ¹⁸O values of the 19 snow pit samples used to assess the models' capability to reproduce observed δ¹⁸O ratios in Greenland are illustrated. In general, COSMO_iso in a 50 km×50 km spatial resolution is able to reproduce the observed isotopic ratios at the snow pit samples and improves the simulation results of ECHAM5-wiso. In both simulations, the δ¹⁸O ratios are high near the coastline and low in central Greenland. But in COSMO_iso_50km, the δ¹⁸O ratios decline more rapidly from the coastline to the inland plateau than in ECHAM5-wiso. The spatial δ¹⁸O differences are consequently more pronounced, and the general overestimation of δ¹⁸O ratios, which occurs in ECHAM5-wiso, is reduced in COSMO_iso_50km. As a consequence, the regional simulation reaches a better agreement with the observations. In COSMO_iso_50km the root-mean-squared error (RMSE) is reduced by 0.98 ‰ over all snow pit samples (the RMSE of ECHAM5-wiso is 2.42 ‰, and the RMSE of COSMO_iso_50km is 1.44 ‰). This RMSE reduction is significant at the 95 % level, assessed by performing a t test. Especially for the snow pit samples for which ECHAM5-wiso exhibits strong deviations from the observed δ¹⁸O values ($\mathrm{1},\mathrm{3},\mathrm{4},\mathrm{6},\mathrm{7},\mathrm{8},\mathrm{9},\mathrm{10},\mathrm{16},\mathrm{17},\mathrm{19}$; see Table 2), a regional downscaling with COSMO_iso_50km reduces the bias considerably (Fig. 4). For these stations, the RMSE of ECHAM5-wiso of 3.09 ‰ is reduced by 1.65 ‰ to 1.44 ‰. But for snow pit samples at which ECHAM5-wiso already has a high agreement with the observations $(\mathrm{2},\mathrm{5},\mathrm{11},\mathrm{13},\mathrm{14})$, COSMO_iso_50km increases the RMSE from 0.34 ‰ to 1.51 ‰.

Figure 3 Yearly mean δ¹⁸O values of COSMO_iso_50km (a) and ECHAM5-wiso (b, interpolated to the COSMO_iso_50km grid) for the period 2008–2014 and the corresponding observations for the 19 snow pit samples (Table 2).

Figure 4 Differences (Δ) between simulated and observed δ¹⁸O values (model minus observation) for the model results of ECHAM5-wiso, COSMO_iso_50km, and COSMO_iso_7km and snow pit samples and top core samples from ice cores from Greenland (averaging periods: COSMO_iso_50km: 2008–2014; COSMO_iso_7km: 2011; observation: 1940–2014). Numbers refer to the different snow pit locations shown in Fig. 3.

Figure 5 shows that these annual biases of the COSMO_iso_50km simulation are not caused by systematic seasonal biases, as for example reported by Sjolte et al. (2011) for RCM simulations in Greenland. Shown are the simulated monthly δ¹⁸O values with COSMO_iso_50km compared to observed monthly δ¹⁸O values in precipitation for the period 2008–2014, collected at Arctic stations of the GNIP data set (Table 1). In general, the modeled δ¹⁸O values for precipitation are in good agreement with the monthly GNIP data. But in the COSMO_iso_50km simulation no systematic overestimation or underestimation of observed isotope ratios is simulated with the RCM. This is true for each season. Neither in winter (low δ¹⁸O values) nor in summer (high δ¹⁸O values) are systematic deviations to the observations simulated. Thus, the seasonal variability in the COSMO_iso_50km results has no systematic impact on the yearly mean δ¹⁸O values and is therefore not the reason for systematic differences between the coarse model results and observations.

Figure 5 Monthly δ¹⁸O simulated with COSMO_iso_50km for the period 2008–2014 versus the corresponding observations for nine GNIP stations (Table 1).

As is visible in Fig. 3, these systematic differences are instead caused by a southward shift of the area of low yearly mean δ¹⁸O values in central northern Greenland in COSMO_iso_50km relative to ECHAM5-wiso. As a result, the simulated δ¹⁸O values in central northern Greenland in COSMO_iso_50km are higher than in ECHAM5-wiso. Since ECHAM5-wiso already has a high agreement with the observed δ¹⁸O values, a model bias is introduced in COSMO_iso_50km, causing the deviations relative to the observations in northern Greenland. But COSMO_iso_50km in general yields an overall improvement in simulating the yearly mean δ¹⁸O values compared to ECHAM5-wiso.

A further downscaling with COSMO_iso to a spatial resolution of 7 km×7 km does not improve the simulation results further. The RMSE reduction with respect to ECHAM5-wiso is 0.83 ‰ (compared to 0.98 ‰ for COSMO_iso_50km). The only exception constitutes the snow pit sample from Renland (1). Here, a considerable model bias in ECHAM5-wiso and COSMO_iso_50km is strongly reduced in COSMO_iso_7km. The coastal area of Renland is characterized by complex terrain and constitutes a special case for isotope-enabled modeling in Greenland. The snow pit sample is located in a transition zone from the homogeneous inland glaciation to the rugged coastline, where the glaciers calve into the sea. Thus, within short distances large differences in altitude and land surface characteristics occur in this region. The isotopic ratios in the snow pit sample are therefore strongly affected by these heterogeneous local conditions, which are insufficiently represented in the coarse model resolution of ECHAM5-wiso. By increasing the spatial resolution with regional climate modeling, the representation of the associated small-scale processes is improved. This leads generally to an improved agreement of the simulation results with observations, as seen for the COSMO_iso_7km run for Renland (Fig. 4).

However, an increase in spatial resolution is also associated with an increased heterogeneity of the surface characteristics and the related small-scale processes, especially in complex terrains. This is because the GCM grid boxes are further divided in smaller RCM grid boxes and consequently both higher and lower values (e.g., altitude) are now included in the respective GCM grid boxes. As a consequence, an additional spatial variability is introduced in the RCM simulations in comparison to the GCM results. Due to uncertainties accompanied by model simulations, this can potentially increase the RCM bias with respect to in situ point measurements, which may actually be closer to the spatially averaged values simulated by the coarse GCM model. This effect can be observed for the SE-Dome ice core (19) in southeastern Greenland. Comparable to the Renland ice core, SE-Dome is located near the coastline. But in contrast to the Renland ice core, an increase in the spatial resolution to 7 km does not further improve the RCM results for SE-Dome. On the contrary, the δ¹⁸O bias is even higher than in the ECHAM5-wiso simulation.

Despite the lack of improvement in the point to grid cell comparison, higher resolved RCM simulations allow the subgrid-scale variability of δ¹⁸O within GCM grid boxes to be simulated and compared to observed δ¹⁸O values. In this way, the inherent uncertainty of in situ measurements, associated with a local micrometeorological variability, can be considered. Thus, in the following sections, snow pit samples are no longer solely compared to the model grid boxes covering the sample locations. Instead, it is investigated whether the δ¹⁸O range of all adjacent RCM grid boxes to a snow pit location is consistent with the observed δ¹⁸O value of the same site. For this, all RCM grid boxes located within the corresponding GCM grid box are included in the comparison with the observations.

3.1.3 δ¹⁸O variability

The spatial isotopic ratio variability of the COSMO_iso_50km grid boxes surrounding the 19 snow pit samples is shown as a box and whisker plot in Fig. 6a. The spatial isotopic ratio variability of the COSMO_iso_7km is shown in Fig. 6b. This spatial isotopic ratio variability includes the δ¹⁸O values of all COSMO_iso (50 km) and COSMO_iso (7 km) grid boxes within the ECHAM5_wiso grid box closest to the snow pit sample. For 14 of the 19 snow pit samples $(\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6},\mathrm{7},\mathrm{8},\mathrm{14},\mathrm{15},\mathrm{16},\mathrm{17},\mathrm{18},\mathrm{19})$ the observed δ¹⁸O values are within the range of the spatial COSMO_iso_50km grid box variability. But for 5 of the 19 snow pit samples (9–13) the spatial isotope range of the COSMO_iso_50km simulation does not fit with the observations. Since these stations are all located in the north of Greenland (Fig. 3), this is likely partially associated with the southward shift of the area of low yearly mean δ¹⁸O values in central northern Greenland in COSMO_iso_50km in comparison to ECHAM5-wiso, as already described in Sect. 3.1.2.

Figure 6 Present-day isotopic ratio variability of the COSMO_iso grid boxes surrounding the 16 snow pit samples for the (a) 50 km and (b) 7 km simulation. The black bar in the box and whisker plot represents the median of the isotope ratio distribution. The box comprises the upper and lower quartile, and the whiskers show the whole distribution. The MPI-ESM-wiso results are shown by the blue dots, and the observed δ¹⁸O values are shown by the red dots.

A downscaling to 7 km does slightly increase the spread of the COSMO_iso results. However, the observed δ¹⁸O values from 5 of 19 snow pit samples still fall outside the range of the modeled COSMO_iso_7km grid box variability (Fig. 6b). Thus, a further downscaling to a spatial resolution of 7 km does not substantially change the simulated isotopic ratio spread within an ECHAM5-wiso grid box. In accordance with the missing benefits of the COSMO_iso_7km simulation and its increased computing time costs, only a COSMO_iso_50km simulation is performed for the mid-Holocene (Sect. 3.2).

The high spatial δ¹⁸O variability in COSMO_iso simulations is also reflected in the spatial δ¹⁸O–temperature slope of the COSMO_iso_50km run (Fig. 7c), a measure that is frequently used to analyze how strongly isotope ratios and surface temperatures are interrelated. The spatial isotope–temperature slope constitutes a linear fit between the simulated δ¹⁸O ratios and the surface temperatures at all COSMO_iso_50km grid boxes within the respective ECHAM5-wiso grid box. The spatial δ¹⁸O–temperature slope is high in central Greenland (1.5–2.0) and moderate at the coastline (∼1.0). In order to better understand these spatial δ¹⁸O–temperature patterns, both quantities affecting the spatial δ¹⁸O–temperature slope, i.e., the spatial δ¹⁸O variability and the spatial temperature variability, are explicitly shown in Fig. 7a and b. These spatial variabilities are calculated as the standard deviation of all COSMO_iso_50km grid boxes within the respective ECHAM5-wiso grid boxes.

Figure 7 Present-day spatial subgrid-scale variability (calculated as standard deviation) of (a) δ¹⁸O and (b) surface temperature, derived from the COSMO_iso_50km grid boxes within the respective ECHAM5-wiso grid boxes for all of Greenland. Present-day (c) spatial and (d) interannual temporal δ¹⁸O–temperature slope for Greenland, based on yearly mean values.

Figure 7a shows that at the coastline, the spatial δ¹⁸O variability of COSMO_iso_50km is considerably increased within the ECHAM5-wiso grid boxes. In central Greenland the spatial isotopic ratio variability is lower. Thus, the simulated spatial δ¹⁸O variability is high in regions where large orographic differences occur within short distances, like the coastal areas of Greenland, and lower for homogeneous terrain like the inland plateau. Nevertheless, widespread areas with higher spatial isotopic variability also occur in the inland plateau of Greenland. This is not the case for the spatial surface temperature variability. In central Greenland, surface temperature variability is very low (Fig. 7b), but near the coastline the spatial surface temperature variability is also high, highlighting how important the land surface characteristics are for the regional temperature variability. The moderate spatial δ¹⁸O–temperature slope at the coastline (Fig. 7c) is therefore a result of a high surface temperature variability in this region counteracting the high δ¹⁸O variability in the slope calculation. There, the correlation between both quantities is consequently high (Fig. S1a in the Supplement, mainly between 0.7 and 0.99). In central Greenland, the spatial δ¹⁸O–temperature slope is further increased due to the relatively high δ¹⁸O variability compared with the surface temperature variability. Therefore, the spatial distribution of δ¹⁸O cannot be solely explained by land surface processes and the associated spatial temperature variability. The additional spatial δ¹⁸O variability consequently must be caused by dynamic atmospheric processes. In this way, isotopic ratios based on atmospheric fractionation processes along the trajectory of an air mass are transported to central Greenland and increase the isotopic variability there.

In order to investigate the temporal interrelations between the isotope ratios and the surface temperature, the interannual temporal δ¹⁸O–temperature slope is calculated for the COSMO_iso_50km simulation, based on the yearly mean δ¹⁸O and surface temperature values (Fig. 7d). The interannual temporal δ¹⁸O–temperature slope is, in contrast to the spatial δ¹⁸O–temperature slope, small all over Greenland, which is in agreement with the results of Sjolte et al. (2011). That means that the interannual δ¹⁸O variability is less pronounced than the interannual surface temperature variability, and thus the correlation between both quantities is also not as strong as for the spatial δ¹⁸O variability and the spatial temperature variability (Fig. S1b). The mean spatial correlation over Greenland is 0.61, compared to 0.25 for the mean temporal correlation. The impact of interannual surface temperature variations on the temporal δ¹⁸O variability in Greenland is therefore not as dominant as for the spatial δ¹⁸O variability.

3.2 Mid-Holocene

3.2.1 Comparison of simulated δ¹⁸O data to ice core data

In contrast to the present-day simulations, for the mid-Holocene COSMO_iso_50km is driven by MPI-ESM-wiso rather than ECHAM5-wiso. While in ECHAM5-wiso oceanic boundary conditions are prescribed by monthly varying sea surface temperatures and sea ice cover, ocean states are calculated internally in the fully coupled atmosphere–ocean Earth system model MPI-ESM-wiso. Systematic deviations between the COSMO_iso_50km simulations for the mid-Holocene and present day, caused by these different forcing approaches, therefore cannot be excluded. For this reason, a comparison of the mid-Holocene δ¹⁸O anomalies to the present-day conditions is omitted, and an analysis is performed for simulated mid-Holocene δ¹⁸O ratios with comparison to observed mid-Holocene δ¹⁸O values.

In Fig. 8a the absolute differences between the simulated MPI-ESM-wiso (blue) and COSMO_iso_50km (green) grid box results and the observed δ¹⁸O ratios at the corresponding ice cores are presented for the mid-Holocene. As in Fig. 6, the spatial isotopic ratio variability of the COSMO_iso_50km grid boxes surrounding the four Greenland ice core samples is shown as a box and whisker plot. MPI-ESM-wiso properly reflects the isotopic ratios of the mid-Holocene from ice core data. For the inland ice cores (GRIP, GISP2, NGRIP), the simulated δ¹⁸O deviates only by about 1 ‰ relative to the observation, and at Renland the deviation is about 3 ‰. For GRIP and GISP2 the MPI-ESM-wiso simulations slightly underestimate the δ¹⁸O ratios, and for NGRIP and Renland the δ¹⁸O values are slightly overestimated.

Figure 8 (a) Mid-Holocene isotopic ratio variability of the COSMO_iso_50km grid boxes surrounding four Greenland ice core samples. In each grid box, the observed ratios derived from the ice cores are subtracted from the simulated δ¹⁸O ratios. The black bar in the box and whisker plot represents the median of the isotope ratio difference distribution. The box comprises the upper and lower quartile, and the whiskers show the whole distribution. The MPI-ESM-wiso (blue dots) and COSMO_iso_50km (green dots) results for the grid points closest to the ice cores are also shown as differences with respect to the observed δ¹⁸O ratios. (b) Anomalies of the MPI-ESM-wiso simulation relative to pre-industrial (PI) conditions, based on an MPI-ESM-wiso PI reference simulation (Cauquoin et al., 2019) are shown as blue points, and the observed anomalies for the mid-Holocene relative to present-day are shown as red points.

For COSMO_iso_50km, the deviation of δ¹⁸O values relative to the observations are opposite in sign compared with MPI-ESM-wiso at all locations except Renland. That means that in GRIP and GISP2, the underestimated δ¹⁸O values in MPI-ESM-wiso are turned into overestimated δ¹⁸O values in COSMO_iso_50km; at NGRIP the overestimation is turned into an underestimation, but the net bias is not reduced. At Renland, the bias is even increased. Thus, by just looking at the absolute biases, the downscaling does not seem to bring an added value to the MPI-ESM-wiso results for mid-Holocene conditions. However, when the spatial isotopic ratio variability within MPI-ESM-wiso grid cells simulated by COSMO_iso_50km isotopic ratios is taken into account, the model results are in agreement with the isotopic ratios of the ice core samples.

In Fig. 8b, the MPI-ESM-wiso model anomalies with reference to the pre-industrial period (PI) conditions, which are based on an MPI-ESM-wiso PI reference simulation performed by Cauquoin et al. (2019), are compared to the measured mid-Holocene-PI δ¹⁸O anomalies of the ice cores. The positive δ¹⁸O anomalies between the mid-Holocene and PI values for both ice core data and MPI-ESM-wiso model results are associated with higher temperatures, especially during the summer, and a reduction in Arctic sea ice during the mid-Holocene (Cauquoin et al., 2019). In Renland and NGRIP simulated anomalies are slightly underestimated, and in GRIP and GISP2 anomalies are slightly overestimated, but overall the biases of the MPI-ESM-wiso mid-Holocene–PI period model anomalies to the observed mid-Holocene–PI period anomalies are small and statistically not significant at the 95 % level for all ice cores (assessed by performing a t test).

3.2.2 δ¹⁸O variability

The fact that, in contrast to the present-day simulations, only four observational data sets are available for the mid-Holocene, makes the assessment of the simulation results difficult. Moreover, with the GRIP and GISP2 ice cores being located very close to each other (Fig. 9), only three local isotope distributions that are clearly different from each other are available. Therefore, in Fig. 9a, the spatial δ¹⁸O variability of the COSMO_iso_50km simulation for the mid-Holocene is illustrated for all of Greenland, which is, in accordance to the analysis of the present-day simulation, again calculated as the standard deviation of all COSMO_iso_50km grid boxes within the respective GCM grid boxes. In general, the δ¹⁸O variability of COSMO_iso_50km in the mid-Holocene is high at the coastline, while it is lower in central Greenland. The Renland ice core is consequently located in an area of a high isotopic variability, and the GRIP and GISP2 ice cores are located in an area of low isotopic variability. But regions with increased isotopic variability also occur in the inland plateau of Greenland. The NGRIP ice core, for instance, is located in such an area of a moderate isotopic ratio variability. The four ice core drill sites are therefore located in three regions of Greenland with substantially different sub-grid isotopic ratio variabilities.

Figure 9 The same as Fig. 7 but for the mid-Holocene. The locations of the ice core samples are shown in green.

The spatial surface temperature variability in the COSMO_iso_50km mid-Holocene simulation is shown in Fig. 9b. The mid-Holocene simulation shows a high spatial surface temperature variability near the coastline and almost no variability in central Greenland. As a consequence, the spatial δ¹⁸O–temperature slope is moderate at the coastline (∼1.0) and high in central Greenland (1.5–2.0; Fig. 9c). Moreover, the interannual δ¹⁸O–temperature slope is very small over Greenland in the mid-Holocene (the mean interannual slope over Greenland is 0.47), although in some regions high temporal slopes are simulated (Fig. 9d). However, the influence of interannual surface temperature variations on the temporal δ¹⁸O variability in the mid-Holocene is, in general, not as dominant as for the spatial variability.

In general, the results of the COSMO_iso_50km mid-Holocene simulation exhibit the same spatial characteristics as for the present-day simulation (Figs. 7 and 9). Comparable spatial patterns are simulated for the surface temperature variability (Figs. 7b and 9b) and the δ¹⁸O variability (Figs. 7a and 9a) within a GCM grid box, although regions of increased δ¹⁸O variability in central Greenland are more widely present in the mid-Holocene run than in the present-day one. In the inland plateau, regions of a high spatial δ¹⁸O–temperature slope are therefore more extensive in the mid-Holocene than under present-day conditions (Figs. 7c and 9c). Nevertheless, the spatial δ¹⁸O–temperature interrelations are comparable in both periods. This is also the case for the temporal variabilities of δ¹⁸O and the surface temperature in the COSMO_iso_50km simulation results for the mid-Holocene and the present day (Figs. 7d and 9d), although both simulations are driven by two different forcing approaches (ECHAM5-wiso nudged to ERA-Interim reanalysis with prescribed monthly varying oceanic boundary conditions vs. the fully-coupled atmosphere–ocean Earth system model MPI-ESM-wiso). This finding indicates that the spatial and interannual δ¹⁸O variability of COSMO_iso_50km within a GCM grid box over Greenland is not strongly dependent on the oceanic boundary conditions.

4 Discussion and conclusions

The results of several global paleo-climate simulations exhibit considerable deviations from the observed regional climate patterns during the Holocene (Braconnot et al., 2012). In the present study, regional climate simulations with an isotope-enabled RCM are, for the first time, performed for Greenland to potentially improve the agreement with climate observations in this region for the mid-Holocene. In a first step, the capability of the isotope-enabled RCM COSMO_iso to reproduce observed isotopic ratios for Greenland is demonstrated.

The COSMO_iso simulation results show that a spatial resolution of 50 km produces reasonable δ¹⁸O values. Especially in regions where the global ECHAM5-wiso model, which has been used to derive necessary forcing fields for the COSMO_iso simulations, deviates strongly from the observed δ¹⁸O values, the RMSE is significantly reduced by 1.65 ‰ for regional climate simulations with COSMO_iso. In complex terrain like the coastal areas of Greenland, the results can be further improved with an additional downscaling to a spatial resolution of 7 km. In such simulations with high spatial resolution, small-scale processes are described in more detail (e.g., Torma et al., 2015; Coppola et al., 2020), and thus the local characteristics at ice core sites are better taken into account (Sturm et al., 2005; Werner et al., 2011). But for the rest of Greenland, highly resolved regional climate simulations do not yield further improvements. For northern Greenland, regional climate simulations with COSMO_iso increase the bias with respect to observations. A comparison of simulated isotope ratios in precipitation with measured values at GNIP stations shows that such deviations between model results and observations are not caused by systematic seasonal biases in the RCM, as it was the case in a similar study by Sjolte et al. (2011) for Greenland. In central northern Greenland, a model bias is instead introduced, due to a southward shift of the area of low yearly mean δ¹⁸O values. All in all, the results of this study show that COSMO_iso is generally able to provide reasonable isotopic ratios for Greenland and that the model can be applied for paleo-climate simulations.

For the mid-Holocene, MPI-ESM-wiso is in good agreement with observed ice core data in Greenland, as already described by Cauquoin et al. (2019). The model bias is, in this context, not further reduced by a downscaling with COSMO_iso. But an increase in the spatial model resolution also leads to an increase in the models' degrees of freedom. This in turn can lead to additional noise and thus a deviating RCM behavior and even an increase in the absolute model bias, as seen for the Renland station.

Another consequence of these increased degrees of freedom in the COSMO_iso simulation is that the spatial variability of the simulated δ¹⁸O ratios is enhanced. This enhanced spatial variability represents the subgrid-scale uncertainty of the driving GCM, which can be derived in a physically consistent way by regional downscaling. By analyzing this subgrid-scale variability, the spatial uncertainties in the comparison between GCM data and point measurements can be considered. In this way, it can be demonstrated that most of the observed δ¹⁸O values lie within the local δ¹⁸O uncertainties of the coarse GCM results. This applies for both the present-day runs and the regional paleo-climate simulations for the mid-Holocene in Greenland. The deviation between the more coarsely resolved GCM results and the more finely resolved observations is therefore potentially caused by the missing representation of important small-scale processes, which are induced by, e.g., the surface conditions or orographic effects over Greenland. Shi et al. (2020), for instance, were able to demonstrate that GCM deficiencies in reproducing the observed water isotope variability in the southeastern Tibetan Plateau are associated with the missing representation of such small-scale processes in coarse GCM simulations.

As δ¹⁸O ratios are used as an indicator for temperatures in past climates (Dansgaard et al., 1969; Masson-Delmotte et al., 2005; Jouzel, 2013), it is important to understand how the presented COSMO_iso simulations might be able to improve these isotope-based temperature reconstructions. In general, the regional surface temperature variability and the regional δ¹⁸O variability show similar patterns for Greenland. In both cases the variability is high at the coast and low on the inland plateau. Similar patterns to those in the mid-Holocene can also be seen for the present-day simulations. These patterns of spatial variability in δ¹⁸O and the surface temperature are in line with the results of Sjolte et al. (2011) for RCM simulations under present-day conditions for Greenland. Based on these patterns of variability, it can be derived that the regional surface temperature variability highly depends on the surface characteristics in Greenland. However, for the regional isotopic ratio variability, this dependence appears to be less pronounced. At the coastline, a clear relationship between surface temperatures and measured δ¹⁸O ratios in ice cores can be deduced, while in central Greenland this relation is weaker. These spatial differences might be explained by the fact that isotope changes are an integrated signal of the mesoscale variability of atmospheric processes (Dansgaard, 1964; Merlivat and Jouzel, 1979; Gat, 1996), which might be partially decoupled from surface temperature changes in homogeneous terrain.

Consistent patterns over Greenland are also modeled for the interannual temporal δ¹⁸O–temperature slope in the mid-Holocene and the present-day simulation. But in comparison to the spatial δ¹⁸O–temperature slope, the interannual temporal δ¹⁸O–temperature slope is rather small. This weaker interannual δ¹⁸O–temperature slope is again in line with the results of Sjolte et al. (2011).

The presented study demonstrates that the isotope-enabled MPI-ESM-wiso–COSMO_iso model chain with realistically implemented stable water isotope fractionation processes constitutes a useful supplement to reconstruct regional paleo-climate conditions during the mid-Holocene in Greenland. By means of such an isotope-enabled GCM-RCM model chain, locally measured isotope ratios in an ice core can be adequately linked to spatially coarse climate model results, and conclusions on the underlying climatic processes leading to these ratios can be drawn in a physically consistent way. This approach might also be very helpful for other isotope-enabled GCMs and understanding their deviations from observed isotope ratios in different paleo-time periods and regions. Particularly in regions in which large differences occur between simulated and observed δ¹⁸O ratios, due to small-scale orographic variations, like parts of Europe and North America (Cauquoin et al., 2019; Comas-Bru et al., 2019), an improved representation of small-scale processes can potentially reduce these biases, and consequently the reconstruction of regional paleo-climate patterns can become more reliable. To test this hypothesis, in follow-up studies more time slices will be simulated with the presented MPI-ESM-wiso–COSMO_iso model chain for different periods and different regions.