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 δ 18 O and observed δ 18 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 δ 18 O ratios from snow pit samples, Global Network of Isotopes in Precipitation (GNIP) stations and ice cores. To our knowledge, this is the ﬁrst 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 (1 .

ment. Despite this lack of improvement in model biases, the study shows that in both periods, observed δ 18 O values at measurement sites constitute isotope ratios which are mainly within the subgrid-scale variability of the global ECHAM5wiso and MPI-ESM-wiso simulation results. The correct δ 18 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 δ 18 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 isotopeenabled 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.
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 cor-responding 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, nearsurface 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.  (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 δ 18 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 isotopeenabled 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.

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 ECHAM5wiso simulation have been already used as boundary conditions for COSMO_iso simulations over Europe by Christner et al. (2018). The simulation outputs from ECHAM5wiso are at a T106 horizontal spatial resolution (1.1 • × 1.1 • ) 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 0.44 • × 0.44 • , 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 0.0625 • × 0.0625 • (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 reanalysisbased 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.
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 cli- matology 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.

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 δ 18 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; Table 2) are calculated as an average of all available δ 18 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). 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 δ 18 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 δ 18 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 alreadymentioned 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 δ 18 O values (blue numbers in Table 2) are calculated as an average of the measured δ 18 O values in ice cores over the period 5.5-6.5 ka.

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  Figure 3 shows the yearly mean δ 18 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 δ 18 O values of the 19 snow pit samples used to assess the models' capability to reproduce observed δ 18 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 ra-tios at the snow pit samples and improves the simulation results of ECHAM5-wiso. In both simulations, the δ 18 O ratios are high near the coastline and low in central Greenland. But in COSMO_iso_50km, the δ 18 O ratios decline more rapidly from the coastline to the inland plateau than in ECHAM5-wiso. The spatial δ 18 O differences are consequently more pronounced, and the general overestimation of δ 18 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 ECHAM5wiso exhibits strong deviations from the observed δ 18 O values (1,3,4,6,7,8,9,10,16,17,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 (2,5,11,13,14), COSMO_iso_50km increases the RMSE from 0.34 ‰ to 1.51 ‰. 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   (Table 1). In general, the modeled δ 18 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 δ 18 O values) nor in summer (high δ 18 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 δ 18 O values and is therefore not the reason for systematic differences between the coarse model results and observations.
As is visible in Fig. 3, these systematic differences are instead caused by a southward shift of the area of low yearly mean δ 18 O values in central northern Greenland in COSMO_iso_50km relative to ECHAM5-wiso. As a result, the simulated δ 18 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 δ 18 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 δ 18 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   (Table 1). 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 δ 18 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 δ 18 O within GCM grid boxes to be simulated and compared to observed δ 18 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 δ 18 O range of all adjacent RCM grid boxes to a snow pit location is consistent with the observed δ 18 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.

δ 18 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 δ 18 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 (1,2,3,4,5,6,7,8,14,15,16,17,18,19) the observed δ 18 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 δ 18 O values in central northern Greenland in COSMO_iso_50km in comparison to ECHAM5-wiso, as already described in Sect. 3.1.2.
A downscaling to 7 km does slightly increase the spread of the COSMO_iso results. However, the observed δ 18 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 δ 18 O variability in COSMO_iso simulations is also reflected in the spatial δ 18 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 δ 18 O ratios and the surface temperatures at all COSMO_iso_50km grid boxes within the respective ECHAM5-wiso grid box. The spatial δ 18 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 δ 18 O-temperature patterns, both quantities affecting the spatial δ 18 O-temperature slope, i.e., the spatial δ 18 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 7a shows that at the coastline, the spatial δ 18 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 δ 18 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 δ 18 O-temperature slope at the coastline (Fig. 7c) is therefore a result of a high surface temperature variability in this region counteracting the high δ 18 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 δ 18 O-temperature slope is further increased due to the relatively high δ 18 O variability compared with the surface temperature variability. Therefore, the spatial distribution of δ 18 O cannot be solely explained by land surface processes and the associated spatial temperature variability. The additional spatial δ 18 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 δ 18 O-temperature slope is calculated for the COSMO_iso_50km simulation, based on the yearly mean δ 18 O and surface temperature values (Fig. 7d). The interannual temporal δ 18 O-temperature slope is, in contrast to the spatial δ 18 O-temperature slope, small all over Greenland, which is in agreement with the results of Sjolte et al. (2011). That means that the interannual δ 18 O variability is less pronounced than the interannual surface temperature variability, Figure 7. Present-day spatial subgrid-scale variability (calculated as standard deviation) of (a) δ 18 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 δ 18 O-temperature slope for Greenland, based on yearly mean values. and thus the correlation between both quantities is also not as strong as for the spatial δ 18 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 δ 18 O variability in Greenland is therefore not as dominant as for the spatial δ 18 O variability.

Mid-Holocene
3.2.1 Comparison of simulated δ 18 O data to ice core data In contrast to the present-day simulations, for the mid-Holocene COSMO_iso_50km is driven by MPI-ESMwiso 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 atmosphereocean 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 δ 18 O anomalies to the present-day conditions is omitted, and an analysis is per-formed for simulated mid-Holocene δ 18 O ratios with comparison to observed mid-Holocene δ 18 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 δ 18 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 δ 18 O deviates only by about 1 ‰ relative to the observation, and at Renland the deviation is about 3 ‰. For GRIP and GISP2 the MPI-ESMwiso simulations slightly underestimate the δ 18 O ratios, and for NGRIP and Renland the δ 18 O values are slightly overestimated.
For COSMO_iso_50km, the deviation of δ 18 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 δ 18 O values in MPI-ESM-wiso are turned into overestimated δ 18 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- 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 δ 18 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 δ 18 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  are shown as blue points, and the observed anomalies for the mid-Holocene relative to present-day are shown as red points.
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 δ 18 O anomalies of the ice cores. The positive δ 18 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 . 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).

δ 18 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 δ 18 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 presentday simulation, again calculated as the standard deviation of all COSMO_iso_50km grid boxes within the respective GCM grid boxes. In general, the δ 18 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.
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 δ 18 O-temperature slope is moderate at the coastline (∼ 1.0) and high in central Greenland (1.5-2.0; Fig. 9c). Moreover, the interannual δ 18 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 δ 18 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 δ 18 O variability (Figs. 7a and 9a) within a GCM grid box, although regions of increased δ 18 O variability in central Greenland are more widely present in the mid-Holocene run than in the presentday one. In the inland plateau, regions of a high spatial δ 18 O-temperature slope are therefore more extensive in the mid-Holocene than under present-day conditions (Figs. 7c and 9c). Nevertheless, the spatial δ 18 O-temperature interrelations are comparable in both periods. This is also the case for the temporal variabilities of δ 18 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 spa-tial and interannual δ 18 O variability of COSMO_iso_50km within a GCM grid box over Greenland is not strongly dependent on the oceanic boundary conditions.

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 δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 O values lie within the local δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 O-temperature slope in the mid-Holocene and the present-day simulation. But in comparison to the spatial δ 18 O-temperature slope, the interannual temporal δ 18 O-temperature slope is rather small. This weaker interannual δ 18 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 paleoclimate 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 isotopeenabled 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 δ 18 O ratios, due to smallscale orographic variations, like parts of Europe and North America 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.