High‐Resolution Nudged Isotope Modeling With ECHAM6‐Wiso: Impacts of Updated Model Physics and ERA5 Reanalysis Data

We present here results of new isotope‐enabled simulations with an enhanced ECHAM6‐wiso model version nudged to the ERA5 reanalyses, at two different spatial resolutions, for the period 1979–2018. The isotopic content of snow on sea ice is considered, yielding surface water vapor with lower isotope ratios over sea ice covered areas, and the kinetic fractionation factors for oceanic evaporation are assumed as independent of wind speed. Also, the supersaturation equation was slightly re‐tuned for a better agreement with the Antarctic isotope observations. In addition to the spatial resolution, the impacts of the improved ECHAM6 model physics and the chosen updated ERA5 reanalyses data set on our simulation results are investigated. For this purpose, detailed comparisons to simulation results obtained from the predecessor ECHAM5‐wiso model nudged to ERA5 reanalyses and from ECHAM6‐wiso nudged to ERA‐Interim are performed. Compared to the ERA‐Interim reanalyses, the nudging to ERA5 does not result in substantial changes of modeled surface isotope values on a global scale but water transport over the tropics clearly changed, with increased precipitation amounts over the Amazonian area and related changes of isotopic contents in water vapor in the high troposphere. An ECHAM6‐wiso simulation at high spatial resolution (0.9° horizontal resolution and 95 vertical levels) generally further improves the modeling of annual mean isotope values. The modeling of temporal isotopic variations at seasonal and (sub‐)daily time scales is also slightly improved for the high spatial resolution configuration. The latter will be very useful for an improved interpretation of various water isotope records.

. They are also used to better understand the mechanisms that control the transport of water vapor (e.g., Bonne et al., 2019;Steen-Larsen et al., 2015, Steen-Larsen, Sveinbjörnsdottir, et al., 2014 in link to weather characteristics like wind directions, relative humidity at the atmosphere-ocean interface or the source of moisture. For several decades, general circulation models (GCMs) equipped with prognostic stable water isotopes greatly contributed to a better understanding of these mechanisms, especially through detailed comparison with isotopic data. Within the last decade, high temporal resolution isotope records in precipitation, water vapor, snow/ice or tree rings have been popularized. At the same time, isotope-enabled GCMs have shown improved performance in terms of (spatial and temporal) resolution and relative computing cost. By these technical developments, isotope model-data comparisons have become more and more important in order to evaluate various processes and contributions determining an isotopic signal recorded in natural archives (e.g., Chiang et al., 2020;Werner et al., 2018) or to quantify the controlling mechanisms on water vapor isotopic composition (e.g., Bonne et al., 2019;Risi et al., 2012;Steen-Larsen et al., 2017). Especially for the modern period, the use of atmospheric GCMs nudged to weather forecast reanalyses has proven to be a powerful way to obtain model outputs under the same weather conditions as the sampling time of the observations. This method has been employed for diverse studies (e.g., Bonne et al., 2019;Butzin et al., 2014;Goursaud et al., 2018;Risi et al., 2010;Servettaz et al., 2020) involving the use of isotope-enabled GCMs nudged to reanalyses data like NCEP, ERA-40 or ERA-Interim (Dee et al., 2011;Kalnay et al., 1996;Uppala et al., 2005).
Here, we present new isotopic simulation results from the latest isotope-enabled version of ECHAM6, ECHAM6-wiso (Cauquoin et al., 2019), nudged to the ERA5 reanalyses over the period 1979-2018. ECHAM6wiso has been updated in several aspects since the last version published in Cauquoin et al. (2019), based on the last findings from water isotope observations campaigns. The purpose of this study is to evaluate how well the new model and reanalyses combination "ECHAM6-wiso nudged to ERA5" can simulate different isotope observations and to separate the potential simulation improvements due to (a) a switch from ECHAM5-wiso (Werner et al., 2011) to ECHAM6-wiso that differs in many physical model aspects (Sections 2 and 3.2), (b) a switch from ERA-Interim to ERA5 and (c) an increase in model resolution from T63L47 (approx. 1.9° × 1.9° horizontal resolution and 47 vertical levels) to T127L95 (approx. 0.9° × 0.9° horizontal resolution and 95 vertical levels).

ECHAM6-Wiso: Background
ECHAM6 (Stevens et al., 2013) is the sixth generation of the atmospheric general circulation model ECHAM, developed at the Max Planck Institute for Meteorology. It consists of a dry spectral-transform dynamical core, a transport model for scalar quantities other than temperature and surface pressure, a suite of physical parameterizations for the representation of diabatic processes, and boundary datasets for externalized parameters (trace gas and aerosol distributions, land surface properties, etc.). ECHAM6 forms the atmospheric component of the fully coupled Earth system model MPI-ESM Mauritsen et al., 2019). The implementation of the water isotopes in ECHAM6 as part of MPI-ESM has been described in detail by Cauquoin et al. (2019), and this model version has been labeled ECHAM6-wiso. Overall, the ECHAM6 model represents the present climate as well as, or better than, its predecessor ECHAM5 (Stevens et al., 2013). Regarding the implemented water isotope module, the ECHAM6-wiso model version used in this study has been updated in several aspects as compared to the initial model release (Cauquoin et al., 2019) to make more consistent the model results with the last findings based on water isotope observations.

Supersaturation Equation
The deuterium-excess (defined as d-excess = δ 2 H − 8 × δ 18 O) in the Antarctic snow is strongly influenced by kinetic fractionation effects occurring during snow formation at very low temperatures (Jouzel & Merlivat, 1984). The supersaturation function S is the main model parameter controlling the strength of these kinetic effects. In our model setup, the empirical supersaturation equation has been slightly changed from S = 1.01−0.0045⋅T cond (Werner et al., 2011) to S = 1.02−0.0045⋅T cond (see Text S1 and Figure S1 in the Supporting Information S1 for details) for a better agreement of ECHAM6-wiso model results with the d-excess/ δ 2 H relationship observed in Antarctic snow (Masson-Delmotte et al., 2008). This update of the supersaturation equation has no significant impact on modeled isotope values at locations other than the poles (see Section 4.2 below). The new formulation of the supersaturation equation in ECHAM6-wiso is still close and comparable to the equations S = 1.0−0.004⋅T cond used in LMDZ-iso (Risi et al., 2010), S = 1.0−0.005⋅T cond used for HadCM3 (Tindall et al., 2009), S = 1.0−0.003⋅T cond parameterized in MIROC5-iso (Okazaki & Yoshimura, 2019) or S = 1.0−0.002⋅T cond defined in iCAM5 (Nusbaumer et al., 2017). Bonne et al. (2019) show the importance of considering the isotopic composition of snow-covered sea ice surfaces to reproduce adequately the isotopic composition in water vapor over sea ice covered areas (see the red and blue dark curves in their Figure 1). Based on their findings, instead of assuming that the isotopic composition of sea ice surface always equals the ocean water just beneath the sea ice (an assumption common in isotope-enabled GCMs), the isotopic composition of sea ice surfaces in ECHAM6-wiso is also reflecting the isotopic composition of snow deposited on this surface (see Equation 6 in Bonne et al., 2019), that is, the sublimation (without fractionation effect) of deposited snow on sea ice influences the isotopic composition of the above surface water vapor. Incorporating this process in ECHAM6-wiso leads to a stronger depletion of surface water vapor over such sea ice covered areas.

Wind Speed Dependency of the Kinetic Fractionation Effect During Oceanic Evaporation
Analyzing two years of continuous surface vapor isotopes measurements made on board of a research vessel, Bonne et al. (2019) did not find any wind speed dependency of the kinetic fractionation during evaporation over the Atlantic Ocean, challenging the standard model approach by Merlivat and Jouzel (1979). Moreover, simulation results of ECHAM5-wiso, the predecessor of ECHAM6-wiso, agree better with these isotope measurements if a constant fractionation factor is assumed. On the other hand, a comparison of ECHAM5-wiso and ECHAM6-wiso results with isotopes in surface vapor data measured at Bermuda, Iceland, NEEM Greenland (Steen-Larsen, Steen-Larsen et al., 2015;Steen-Larsen, Sveinbjörnsdottir, et al., 2014) and from cruises over the Atlantic (Benetti et al., 2017) yields slightly different results. It can be partly explained by the fact that the data stems from different isotope analyzers applied during different cruises (with different calibration routines). For this study, we assume a constant kinetic fractionation factor with mean values in-between rough and smooth wind regime ( ; Benetti et al., 2014) instead of the standard model approach by Merlivat and Jouzel (1979). We evaluate this approach with the isotope observations in precipitation in Section 4.2 and in Supporting Information S1. This approach is consistent with the main conclusion of Bonne et al. (2019). So, in this study we use the above-mentioned constant coefficient between rough and smooth wind regime suggested by Benetti et al. (2014) for our nudged ECHAM6-wiso simulations. Uncertainties remain for the choice of this parameter and additional experimental data sets are required for any further improvement of the chosen value.

ERA Reanalyses
The ERA-Interim reanalyses data, whose production was discontinued at the end of August 2019, has been superseded by the follow-up hourly reanalysis data set ERA5 (Hersbach et al., 2020). The data processing is carried out by European Center for Medium-Range Weather Forecasts (ECMWF), using ECMWF's Earth System model IFS (Integrated Forecasting System), cycle 41r2. The reanalyses are currently available for the period 1979-2021, and a preliminary version of "back extension" data covering the period 1950-1978 has been recently released. Compared to ERA-Interim, the temporal and spatial resolution has been improved: from 6-hourly in ERA-Interim to hourly in ERA5, and from 79 km in the horizontal dimension and 60 levels in the vertical to 31 km and 137 levels in ERA5. In addition to the resolution, the IFS model used for producing the ERA5 reanalyses data has a much improved troposphere representation, a better global balance of precipitation and evaporation, a better representation of precipitation over land in the deep tropics and more consistent sea surface temperature (SST) and sea ice (SIC) data (Hersbach et al., 2020). In addition of the sub-daily reanalyses data used for nudging, monthly mean ERA-Interim (Dee et al., 2011) and ERA5 (Hersbach et al., 2019) data are also employed for analyses shown in Section 4.4.

Simulation Setups
All ECHAM5/6-wiso simulations performed for this study have their 3D-fields of temperature, vorticity and divergence as well as their surface pressure field nudged toward the chosen ERA data every 6 hr for the period 1979-2018. The humidity field is not nudged in our simulations. A spin-up by repeating consecutively 30 times the weather conditions of the year 1979 (i.e., 30 simulation years) has been performed at the start of all the simulations. The model outputs are extracted every 6 hr, allowing model-data comparisons at high-temporal resolution. The orbital parameters and greenhouse gases concentrations have been set to the values of the corresponding model year. The monthly mean SST and sea-ice fields from the corresponding ERA reanalyses have been applied as ocean surface boundary conditions for our different experiments, as well as a mean δ 18 O of surface seawater reconstruction from the global gridded data set of LeGrande and Schmidt (2006). As no equivalent data set of the δ 2 H composition of seawater exists, the deuterium isotopic composition of the seawater in any grid cell has been set equal to the related δ 18 O composition, multiplied by a factor of 8, in accordance with the observed relation for meteoric water on a global scale (Craig, 1961). This assumption implies a surface seawater d-excess of 0‰.
In order to evaluate the effects of the model spatial resolution, two ECHAM6-wiso simulations at T63L47 (approx. 1.9° × 1.9° horizontal resolution and 47 vertical levels) and T127L95 (approx. 0.9° × 0.9° horizontal resolution and 95 vertical levels) spatial resolutions have been performed (Table 1). To evaluate the impacts of ERA-Interim versus ERA5 reanalysis data set for the nudging, a simulation with ECHAM6-wiso nudged to ERA-Interim has been performed, too. The change of model version is also investigated by using the predecessor ECHAM5-wiso instead of ECHAM6-wiso, run at T63L31 spatial resolution. Many physical processes (e.g., radiation scheme, albedo, cloud cover, soil water scheme) as well as the vertical discretization of the atmosphere have been largely updated from the version 5 to version 6 of ECHAM (Stevens et al., 2013). Because the production of the ERA-Interim data set has been stopped in August 2019, all the analyses are performed for the period from January 1979 to December 2018.

Observational Data
To evaluate our simulations, we use different isotopic datasets that we briefly describe here. For the isotopic content of precipitation at global spatial scale, we use the Global Network of Isotopes in Precipitation (GNIP) database, available through the International Atomic Energy Agency (IAEA, 2020). For the annual mean values, we use the same subset of 70 stations from this database as in Cauquoin et al. (2019). In addition to the GNIP data, we use here the same selection of Greenland and Antarctic ice cores (five and 10 locations, respectively) as in Cauquoin et al. (2019) to compare the measured δ 18 O values of polar precipitation with our model results (see Figure 1 for a map of all selected positions). For the seasonal variations in water isotopes, temperature and precipitation amount, we selected six GNIP stations with sufficient time coverage within the 1979-2018 period and representative of different climate conditions (e.g., temperate, tropical, polar, or continental climate): Vienna (central Europe, spanned period: 1979, n = 464), Valentia (coastal Europe, 1979, n = 426), Reykjavik (Arctic, 1992, n = 282), Belem (Amazonian tropical area, 1979-1987, n = 83), Halley Bay (Antarctic coast, 1979 and Ankara (Eastern Mediterranean, 1979. Due to the high temporal resolution of our simulations outputs, we can also compare the model results to several records of (sub-)daily isotopic variations in surface water vapor for different climatic zones: Ankara, Turkey (GNIP); Niwot Ridge, USA (Berkelhammer et al., 2016); Mase, Japan (Wei et al., 2015(Wei et al., , 2016; and Summit, Greenland (Bailey et al., 2015). For Summit data, we selected the data at the height of 19 m. The other stations Ankara, Mase and Niwot Ridge had inlet heights of 5, 2 and 24.3 m above the ground, respectively. The isotope data series with a temporal resolution higher than 6 hr, like Mase, Summit and Niwot Ridge records, have been smoothed and resampled to the ECHAM5/6-wiso temporal resolution (i.e., 6 hr) for the analyses. The comparisons to the observations are made by extracting the modeled values at the nearest grid cell of the simulation if not stated otherwise. The model outputs related to surface water vapor are extracted from the first vertical model layer above the surface, which is approximately 60 m thick. Finally, the ERA5 and ERA-Interim reanalysis data sets are also used to evaluate the changes in modeled water vapor, its transport and its isotopic content in our different simulations.

Simulated Annual Mean δ 18 O p and Temperature Values
We evaluate here the global distribution of the δ 18 O p in precipitation (δ 18 O p ) and its relationship with near surface air temperature from our different simulations. Figure 1 shows the comparison of the modeled δ 18 O p with the isotope observations from the GNIP data set and the added ice cores compilation ( Figure 1a). A very good agreement with the observations is found for all the simulations (δ 18 O p model-data slope of 0.85 at least for all the simulations, Table 2), with a slight improvement for the ECHAM6-wiso simulation results The root mean square error (RMSE) is also improved in the same way. The best model-data agreement is found with the higher resolution E6_HR_ERA5 simulation, shown by the smaller RMSE (2.1‰) as compared to the other three simulations (2.4-2.8‰, see Table 2). The δ 18 O p -temperature relationships (for annual mean surface temperatures below 20°C) retrieved from our simulation results are in very good agreement with the observations with modeled gradients of 0.62-0.63‰°C −1 (r 2 = 0.97), very close to the observed one (0.66‰°C −1 , r 2 = 0.95, Table 2). There is no significant difference between the modeled δ 18 O p -temperature gradients of the four simulations. After having checked that our simulations reproduce adequately the observed δ 18 O p distribution, we focus now on the differences in modeled δ 18 O p between the four ECHAM5/6wiso simulations. The differences in modeled δ 18 O p between ECHAM6-wiso and ECHAM5-wiso ( Figure 1d) are large, especially over sea ice covered areas with anomalies reaching more than 4‰. This is a result of the isotopic composition of snow deposited on sea ice surfaces that is taken into account in ECHAM6-wiso (Section 2.3). Higher δ 18 O p values over land (between 1 and 2‰) in ECHAM5-wiso compared to ECHAM6-wiso can be explained by numerous changes in land surface processes, including the simulation of the soil hydrology within five soil layers instead of one (Hagemann & Stacke, 2015). Differences in temperature or precipitation rate can explain the δ 18 O p anomalies over the poles and tropical oceans, respectively. They can be due to substantial updates in radiation scheme, cloud cover and albedo (Stevens et al., 2013).
Finally, the change of spatial resolution ( Figure 1e) uniformly affects the modeled δ 18 O p over ocean covered areas (decrease of down to −1‰ when the ECHAM6-wiso spatial resolution is increased). Over land areas, δ 18 O p is increased by more than 1‰ in the East Asian monsoon area and South America. Stevens et al. (2013) suggested that the representation of tropical cloudiness is improved at higher resolution. Lower δ 18 O p values are simulated in E6_HR_ERA5 compared to E6_LR_ERA5 over elevated areas like the Himalayas or the Andes due to the improved topography in E6_HR_ERA5.

Seasonal Cycles of T, P, δ 18 O p , and d-Excess
The simulations of the seasonal variability of δ 18 O p , d-excess in precipitation, the 2m air temperature and precipitation are evaluated for six GNIP locations (Section 3.3). The model-data Pearson's correlation coefficients (r) of monthly mean records are shown in Figure 2, and their values as well as the corresponding  linear correlation gradients are reported in Table S2. Except for the Belem GNIP station, the timing and the amplitude of seasonal temperature variations are simulated consistently with the observations ( Figure S3 in the Supporting Information S1) as shown by a high r value (Figure 2 and Table S2, r ≥ 0.95), which is expected considering that the 3D-field of temperature is nudged in these simulations. At Belem, in addition to the lower model-data correlation for temperature by a factor of more than 2 compared to the other GNIP stations, ECHAM5/6-wiso simulates larger variations of temperature at Belem (colored curves in Figure S3 in the Supporting Information S1, standard deviation σ of at least 1.5°C) compared to the observations (σ = 0.55°C). Belem station is strongly influenced by surrounding tropical sea surface condition. An underestimation of this marine influence in ECHAM5/6-wiso could explain the model biases at Belem location (see Section 4.4). The seasonal variability of δ 18 O p in precipitation is correctly simulated for all the GNIP stations, expect for Belem ( Figure S3 in the Supporting Information S1), with r > 0.6 ( Figure 2). Despite a correct representation δ 18 O p seasonal variability at Halley Bay (Antarctic coast), austral summer δ 18 O p values are underestimated by 5‰ ( Figure S3 in the Supporting Information S1). The model-data correlation in d-excess monthly variability is over 0.4 for Reykjavik (Iceland), Valentia and Ankara, but strong deviations to the observations exist for the other locations. Except for Belem (0.11 ≤ r ≤ 0.57) and Halley Bay (0.26 ≤ r ≤ 0.29), the seasonal timing of precipitation is globally well reproduced at numerous stations (0.62 ≤ r ≤ 0.84) but with biased seasonal amplitude of modeled precipitation at Ankara and Vienna (∼30 mm month −1 in observations against more than 80 mm month −1 in simulations, Figure S3 in the Supporting Information S1). We conclude that ECHAM5/6-wiso reproduces correctly the seasonal variability of near surface air temperature, precipitation and the related water isotope signals. The impacts of reanalyses used for the nudging, model version and resolution on the modeled seasonal variability are now investigated.
For all the stations and all the considered variables, differences between the model-data correlation coefficient values of E6_LR_ERAI and E6_LR_ERA5 simulations (red and blue markers in Figure 2 for ERA-Interim and ERA5, respectively) are less than 0.1 (Table S2). This statement includes Belem too, a tropical  Table 1 for the simulation names): E6_LR_ERAI (red), E5_LR_ERA5 (green), E6_LR_ERA5 (blue) and E6_HR_ERA5 (orange).
area where ERA5 reanalyses are improved compared to ERA-Interim products. So, the modeled seasonal variabilities are influenced very little by the two different ERA reanalyses data sets used for the nudging. We investigate now the effects of switching from ECHAM5-wiso to ECHAM6-wiso on our results (green and blue markers in Figure 2, respectively). For the δ 18 O p signal in precipitation, the agreement with the observations is improved by switching from ECHAM5-wiso to ECHAM6-wiso for the polar GNIP stations Reykjavik (r = 0.66 and 0.69, respectively) and Halley Bay (r = 0.43 and 0.62, respectively) that are both influenced by nearby oceanic water sources. On the other hand, some deteriorations of model-data agreement for the GNIP stations Valentia, Ankara and Belem are also noticed. Concerning the d-excess, a better agreement with the observations is found for Ankara (r = 0.39 and 0.44 for E5_LR_ERA5 and E6_LR_ERA5, respectively), while a worse one is obtained for Vienna (r = 0.33 and 0.20, respectively). At Belem, the seasonal amplitude of modeled d-excess is lower in ECHAM6-wiso than in ECHAM5-wiso (standard deviations σ = 2.3 and 4.0‰ respectively, blue and green curves in Figure S3 in the Supporting Information S1). The modeled mean d-excess value is higher in ECHAM6-wiso than in ECHAM5-wiso (7.14 and 4.76‰, respectively), in better agreement with the observed mean value (8.78‰), though. Smaller seasonal variations of modeled 2m air temperature in ECHAM6-wiso than in ECHAM5-wiso (σ = 1.5 and 2.2°C, respectively) are found too, which are also in better agreement with the observations (σ = 0.55°C). The lower modeled mean temperature in the E6_LR_ERA5 simulation than in E5_LR_ERA5 (25.6 and 27.4°C, respectively) at Belem is also more consistent with the observed one (T 2m = 26.0°C). From all these results, we conclude that the change of ECHAM5/6-wiso version yields a mixed set of results in terms of model-data agreements.
As a next step, the increase of spatial resolution, in both the horizontal and vertical dimension, is investigated (orange markers in Figure 2). Compared to the E6_LR_ERA5 simulation (blue markers in Figure 2), the E6_HR_ERA5 simulation results in clear improvements in the model-data agreement for δ 18 O p at the stations Vienna, Valentia and Reykjavik (r = 0.84, 0.73 and 0.74 for respectively, against r = 0.82, 0.67 and 0.69 for E6_LR_ERA5 simulation). This is probably related to the improved model-data agreement in terms of temperature and precipitation variations (blue and orange markers in Figure 2). The simulation of d-excess variations in the E6_HR_ERA5 simulation is in better agreement with the observations at Belem (r = 0.37 against r = 0.21 in E6_LR_ERA5). At Ankara, the correlation coefficient is slightly deteriorated (r = 0.44 and 0.41 in E6_LR_ERA5 and E6_HR_ERA5, respectively) but the model-data linear correlation gradient is significantly improved (from 0.33 to 0.45, Table S2). For the other stations like Vienna, Halley Bay and Reykjavik, the model-data agreement of the E6_HR_ERA5 simulation with the isotopic observations is comparable to the one obtained with E6_LR_ERA5. We deduce from this analysis that the higher spatial resolution T127L95 globally improves the seasonal variability of both standard variables (i.e., temperature and precipitation) and isotopic contents of precipitation.
It is also worth to note that the introduced update of the supersaturation equation in this ECHAM6-wiso model version (see details in Section 2.2, above) has no significant impact on modeled isotope values at locations other than the poles. As seen in Figure 3, the modified supersaturation equation in ECHAM6-wiso changed the modeled δ 18 O and d-excess in precipitation values in Halley Bay (Antarctic coast) by +0.2 and −0.9‰, respectively (blue and red lines in Figure 3). For all other analyzed locations, the changes in δ 18 O and d-excess are negligible. Figure 3 also displays the agreement of the chosen constant kinetic fractionation factor during the evaporation over the ocean with the previous wind-speed dependent fractionation factor suggested by Merlivat and Jouzel (1979) (purple and red lines, respectively. Statistics reported in Table S1 in the Supporting Information S1). These results show that using no wind speed dependency for the definition of kinetic fractionation factors during the evaporation over the ocean does not deteriorate the modeling of the water isotopes. In contrast, a constant factor valid for a smooth wind regime leads to larger disagreement of d-excess values at all selected GNIP stations except Valentia (green line).

Isotopic Content of Surface Water Vapor
We analyze in this section the variations of modeled  Figure S4 in the Supporting Information S1). The model-data correlation coefficient values of the surface specific humidity variations are at minimum 0.92, 0.96 and 0.66 for Summit, Mase and Niwot Ridge station, respectively.
We analyze here the impacts of our simulation setups on modeled δ 18 O v in surface water vapor (Figures 4 and S5 in the Supporting Information S1). The model-data correlation coefficients range between 0.48 (E6_ LR_ERA5 and E6_LR_ERAI simulations at Niwot Ridge) and 0.76 (E6_HR_ERA5 simulation at Summit) at the considered stations. We notice some general model-data biases like overestimations of modeled δ 18 O v at Ankara, Mase and Niwot Ridge during winter time (DJF period in Figure 4) by 5, 4 and 7‰ respectively. The E6_LR_ERAI and E6_LR_ERA5 simulation results have very close model-data agreements (red and blue lines in Figure 4, respectively), as shown by the differences in model-data correlation coefficient values lower than 0.02 (Table S3). No superior simulation can be detected for the ECHAM5-wiso versus ECHAM6wiso comparison (green and blue lines in Figure 4, respectively), neither. The model-data agreement is better with ECHAM6-wiso at Summit because model results from ECHAM5-wiso spread too much (first row in Figure 4, RMSE = 7.57 and 9.82‰ for E6_LR_ERA5 and E5_LR_ERA5, respectively). On the other hand, ECHAM5-wiso is better at reproducing water isotope variations at Niwot Ridge than ECHAM6-wiso (r = 0.48 and 0.63, respectively). For Ankara and Mase, the ECHAM5-wiso and ECHAM6-wiso simulation results are not distinguishable (Ankara and Mase). The increase of ECHAM6-wiso spatial resolution improves slightly the δ 18 O v model-data agreement at Summit and Niwot Ridge (orange lines in Figure 4, r = 0.76 and 0.54 respectively against r = 0.71 and 0.48 for E6_LR_ERA5 simulation) due to a small bias reduction in June to November months, but deteriorates it at Ankara (r = 0.63 and 0.72 for E6_HR_ERA5 and E6_LR_ERA5, respectively) and Mase (r = 0.55 and 0.62, respectively). We conclude that the change of reanalyses for the nudging, the switching of model version or the increase in spatial resolution does not improve clearly the modeled daily variations of δ 18 O v .
Next, we investigate the impacts on d-excess v in surface water vapor. The d-excess v values modeled by ECHAM5/6-wiso have a negative bias for winter months at Ankara (12‰ in DJF and 5‰ in MAM) and Mase (10 and 9‰ in DJF and MAM, respectively) ( Figure 5). Similarly to Nusbaumer et al. (2017), ECHAM5/6wiso overestimates d-excess v values at Niwot Ridge especially during summer with a model bias reaching more than 10‰ in JJA. We obtain a better agreement with the observations for three of the 4 sites with the higher spatial resolution simulation E6_HR_ERA5 (r = 0.31, 0.21 and 0.12 for Summit, Ankara and Niwot Ridge sites, respectively, see Table S3), showing that higher spatial resolution is an important aspect for a better modeling of d-excess v variations. Switching from ERA-Interim to ERA5 for the nudging does not improve significantly the agreement with the d-excess v observations, regardless of the site considered as shown by the very similar KDE signals in Figure 5 and the differences of less than 0.01 in respective model-data correlation coefficient values. On the other hand, we obtain different model-data agreements according to the used ECHAM version (green and blue lines in Figure 5). Improvements are seen for Summit and Mase sites due to less model spread in boreal winter and a reduction of d-excess v bias (3‰ in JJA), respectively. Deteriorations for Niwot Ridge are noticed especially during DJF with an increase of model bias by 2‰. For the Mase site, the model-data correlation coefficients are close to zero for ECHAM5-wiso and E6_HR_ERA5 simulation while they are 0.26 for both ECHAM6-wiso LR setups. The fact that the nearest grid point of this site is an ocean cell in these ECHAM6-wiso LR setups while this same cell is mainly land in the other simulations (52% and 100% land for the ECHAM5-wiso and E6_HR_ERA5 simulations, respectively) could explain such differences. As an additional test on the impact of higher resolution on our model results, we make the average of the 25 grid cells centered on each station location instead of taking the nearest grid point only. With this method, the model-data correlations in d-excess v are improved for both E6_LR_ERA5 and E6_HR_ERA5 at Summit (r = 0.34 and 0.41, respectively) and Ankara (r = 0.18 and 0.26, respectively) and do not change significantly at Niwot Ridge. For Mase station, the model-data correlation is deteriorated with E6_LR_ERA5 simulation (from 0.26 to 0.21) while it is improved with E6_HR_ERA5 (from −0.09 to 0.19), confirming the importance of the land-sea mask. According to all these results, we conclude that higher resolution is an important factor to improve the representation of the daily variations of d-excess v in surface water vapor. However, we note that this conclusion is based on the comparison of the selected four locations with highly resolved temporal isotope variations in vapor, only.
Next, we compare the global changes in modeled d-excess between our different simulations. Figure 6a shows the modeled d-excess p of precipitation in E6_LR_ERA5 simulation. The d-excess p values range between 0 and 20‰. The largest values are found in the northern part of the Sahara and in a 25-45°N band going from Saudi Arabia to the Himalaya. The lowest values are in dry areas like India and southern Sahara, and over the Southern Ocean (between 2 and 6‰), which is an area with large net freshwater input. For the Antarctic continent, the observed contrast between the low values of d-excess p in the west and high values in the east (Cauquoin et al., 2019) is captured by the model (less than 5‰ and between 8 and 11‰, respectively). The choice of the reanalyses data set for nudging (Figure 6b) impacts the global d-excess p values only slightly. Main changes are in the tropics with higher d-excess p values by 1‰ over western Brazil and Africa around the equator when switching to ERA-Interim instead of ERA5. In contrast, lower values by 0.5‰ are located at tropical latitudes over the Atlantic and Eastern Pacific Oceans. These impacts on modeled d-excess p values are caused by change in water vapor transport in the tropics, investigated in more detail in the following Section 4.4. Moreover, modeled d-excessp is lower in E6_LR_ERAI compared to E6_LR_ERA5 by more than 1‰ over northern Sahara. The switch from ECHAM version 5 to version 6 ( Figure 6c) results in a change of the modeled d-excess p values by 4‰ at maximum. Lower values in E5_LR_ERA5 than in E6_LR_ERA5 over sea ice covered areas (differences between −1 and −4‰) are due to the update for the snow deposited on sea ice surfaces (see Section 2.3). Higher values of d-excess p in Antarctica could be related to changes in temperature, precipitation and the changed supersaturation function between ECHAM5wiso and ECHAM6-wiso. Land surface processes have a clear impact on d-excess p values with lower values (between −1 and −2‰) in Eurasia and America when using ECHAM5 instead of ECHAM6. Higher values by more than 4‰ in E5_LR_ERA5 over Northern Sahara and Saudi Arabia are noticed, too. Values of d-excess p are affected over the Atlantic Ocean and Eastern Pacific Ocean because of the changes in precipitation, caused by diverse updates between ECHAM5 and ECHAM6 (see Section 3.2). The increase of spatial resolution (Figure 6d) reduces the d-excess p by 1.5‰ in Antarctica. A better representation of the topography could explain these differences. Higher values of d-excess p by up to 2‰ are located over Australia and lower values by more than 2‰ are found over southern Sahara. Finally higher d-excess p values by around 0.75‰ in E6_HR_ERA5 compared to E5_LR_ERA5 are located over ocean covered areas at ±25° of latitude around the equator, where the freshwater export is large. To conclude, despite the differences in modeled d-excess p values between our simulations, the agreements with the observations are close to each other for all simulations (model-data slope equal to 0.41 for E6_LR_ERA5 and E6_LR_ERAI, 0.42 for E5_LR_ERA5 and 0.43 for E6_HR_ERA5).

ERA-Interim Versus ERA5 Nudging Fields
The model-data comparisons for surface variables at the selected GNIP and other station sites do not show a clear overall ECHAM6-wiso model result improvement when switching from ERA-Interim to ERA5 for the nudging. However, by analyzing all available simulation results on a global scale instead of the selected GNIP locations, we can identify some regions with stronger changes in isotope-controlling variables, like temperature and precipitation amount, between ERA-Interim and ERA5. In addition, substantial differences in the spatial distribution and transport of water vapor in the atmosphere can be detected for these two reanalyses data sets.
As shown in Figure 7, 2m temperatures, precipitation rate and δ 18 O p modeled by ECHAM6-wiso change only little on a global scale when using ERA-interim or ERA5 as the prescribed nudging fields. Global mean modeled differences for δ 18 O p , 2m air temperature and precipitation rate are −0.21‰, −0.04°C and +0.95 × 10 −3 mm day −1 between E6_LR_ERA5 and E6_LR_ERAI, respectively. The temperatures over the ocean areas are very similar for the two simulations, with the only exceptions in the polar regions where absolute anomalies of temperature can reach 1.5°C (E6 _LR_ERA5 minus E6_LR_ERAI). The 2m air temperatures over land areas are globally lower when using ERA5 reanalyses for nudging. Negative temperature anomalies of more than 1.2°C are simulated over the northern Amazonian area (Figure 7b). On the contrary, positive temperature anomalies of more than 1°C around 10°N are detected in Africa and South America. For precipitation, strong positive E6 _LR_ERA5 -E6_LR_ERAI anomalies, varying between 0.7 and 2.3 mm day −1 , are found in Northern Amazonia (Figure 7c). On a global scale, the rainfall is substantially different over the 15°N-15°S latitudinal band with anomalies between −1.2 and +1.3 mm day −1 , related to a different position of the InterTropical Convergence Zone (ITCZ) between ERA5 and ERA-Interim (Hersbach et al., 2020). These differences in temperature and precipitation have impacts on the modeled δ 18 O p values (Figure 7a). Higher modeled temperatures in E6_LR_ERA5 compared to E6_LR_ERAI over the Arctic Ocean above 80°N (Figure 7b) lead to positive anomalies in δ 18 O p values reaching +1.6‰ (Figure 7a). Similarly, negative δ 18 O p anomalies are found over Eurasia and Greenland (down to −0.8‰, Figure 7a) due to lower mean land temperatures in E6_LR_ERA5 compared E6_LR_ERAI (Figure 7b). Concerning the precipitation change effect on modeled δ 18 O p values, positive δ 18 O p anomalies are located over a tropical latitudinal band (up to +0.9‰, Figure 7a) in link to lower modeled precipitation over this area (Figure 7c). Interestingly, the positive δ 18 O p anomalies over eastern Brazil do not seem to be related to higher temperature or lower precipitation, but could be related to change in water vapor transport (see below). In contrast to the positive δ 18 O p anomalies over eastern Brazil, negative δ 18 O p anomalies are simulated over the western Amazonian area (between −0.2 and −0.8‰, Figure 7a) and in the south of the continent, related to the enhanced modeled precipitation and lower temperatures (Figures 7b and 7c).
One possible process causing the modeled δ 18 O p anomalies in the tropics is a different large-scale water vapor transport in ECHAM6-wiso when ERA5 or ERA-interim is used for nudging. Thus, we compare the amounts and transport directions of vertically integrated water vapor both from the ERA reanalyses data directly and from the related ECHAM6-wiso model results (Figure 8). Because the ECHAM6-wiso simulations are nudged to ERA reanalyses, similar modeled wind patterns and thus similar modeled water vapor transport patterns are expected. The water vapor transport over the tropical Atlantic, from the Gulf of Guinea to Brazil, is enhanced in ERA5 compared to ERA-Interim (Figures 8a and 8b) by around 15 kg m −1 s −1 . Once arrived at the western side of the Amazonian area, part of this water vapor is then deflected to the south along the Andes. Because this pattern is stronger in ERA5 than in ERA-Interim, the water vapor amount over the Amazonian rainforest is higher in ERA5 by ∼2 kg m −2 (Figure 8b). Similar results are found when comparing E6_LR_ERA5 and E6_LR_ERAI simulations instead of the reanalyses (Figures 8c and 8d). We conclude that the stronger transport of water vapor enriched in heavy isotopes (from the marine influence) to Brazil due to nudging to ERA5 instead of ERA-Interim can explain the higher modeled δ 18 O p values over the eastern Amazonian area (Figure 7a) in the E6_LR_ERA5 simulation. Also, we note that the concentration of modeled vertically integrated water vapor over Amazonia in E6_LR_ERA5 is lower than in ERA5 reanalyses (Figures 8a and 8c) by ∼3 kg m −2 . This underestimation by ECHAM6-wiso could explain partly the model deficit in simulating precipitation at Belem station (Section 4.2). The increase in ECHAM6-wiso modeled concentration of water vapor over Amazonia by the change in nudging from ERA-Interim to ERA5 is equal to ∼1 kg m −2 (Figure 8d), two-times lower than the increase found between ERA5 and ERA-Interim reanalyses, directly. This smaller change might be caused by different formulations of atmospheric hydrological processes (e.g., evaporation, cloud formation, condensation) between the IFS model used for producing the ERA renalyses data set and the ECHAM5/6-wiso model. The ERA5 data also indicates an enhanced water vapor transport in the tropical Indo-Pacific region as compared to ERA-Interim (+20 kg m −1 s −1 ). However, there exists no difference in the vertically integrated water vapor content between both reanalyses at this location. For the ECHAM6-wiso model, the water vapor amount over the Northern tropical Indo-Pacific in the E6_LR_ERA5 simulation is slightly lower than in the E6_LR_ERAI simulation by 0.75 kg m −2 (coincident with negative modeled precipitation anomalies, see Figure 7c). However, both simulations underestimate the total amount of water vapor over this area as compared to the reanalyses result (e.g., by 5 kg m −2 in E6_LR_ERA5 simulation compared to ERA5; Figures 8c and 8a, respectively). We conclude from this analysis that clear differences exist between the reanalyses and our ECHAM6-wiso simulations, in terms of vertically integrated water vapor amount especially in the Northern tropical Indo-Pacific area and Amazonia (drier in ECHAM5/6-wiso than in ERA reanalyses), and in terms of transport strength over the tropical oceanic areas. The representation of related hydrological processes like the formation and distribution of clouds or the convective parameterizations could partly explain these discrepancies (Mauritsen  Stevens et al., 2013), with potentially resulting biases in the spatial distribution of the water isotope values. For an improved understanding of these discrepancies, we next consider the vertical distribution of water vapor and of its isotopic content in both the reanalyses and ECHAM6-wiso simulations. As largest differences are found in the tropics, we restrict the analyses on meridional mean vertical tropical profiles of specific humidity over a latitudinal band between 15°S and 15°N ( Figure 9).
The differences in specific humidity between ERA5 and ERA-Interim (Figure 9a) are strong and homogenous along the pressure levels, with absolute values of specific humidity anomalies often over 10 −4 kg kg −1 . Compared to ERA-Interim, the tropical specific humidity in ERA5 is higher at the surface, except over Amazonia (70°W) and Africa (20°E), lower around 870 hPa, higher between 750 and 600 hPa and slightly lower in the high troposphere (between −0.5 and −1 × 10 −4 kg kg −1 ). The strongest ERA5 -ERA-Interim anomalies (more than 5 × 10 −4 kg kg −1 ) are over Amazonia and Africa, i.e., over tropical land areas where ERA5 has been much improved (Hersbach et al., 2020). ERA5 and ERA-Interim are based on the IFS model cycles Cy41r2 and Cy31r2, respectively. ERA5 thus benefits from a decade of developments in model physics, including upgrades of the large-scale cloud and precipitation scheme, and changes in the parameterization of convection (Hersbach et al., 2020). The anomalies in modeled specific humidity between E6_LR_ ERA5 and E6_LR_ERAI simulations (Figure 9b) are not similar to the ones between the ERA reanalyses in several aspects. This can be explained by the already mentioned differences in simulated hydrological processes as well as the fact that the modeled ERA humidity fields are not nudged for these ECHAM5/6wiso simulations. Negative water vapor anomalies are located around the 800 hPa pressure level, except at 70°W (over Amazonia). Negative anomalies are found over Africa (20°E) and Indo-Pacific area (more than 2 × 10 −4 kg kg −1 ), but no positive anomalies higher than 2 × 10 −4 kg kg −1 are found in lower pressure levels (i.e., higher in altitude) contrary to the detected specific humidity anomalies between ERA5 and ERA-Interim. Positive anomalies in modeled water amounts are located over the Amazonian area at around 500 hPa like in the analyzed reanalyses anomalies, but with the anomalies in modeled water amounts being more than three-times lower than the reanalyses anomalies. Negative anomalies are found over the Pacific Ocean around 90°W, and only slightly negative anomalies in modeled specific humidity lower than 10 −4 kg kg −1 Figure 9. Tropical vertical profiles (meridional mean between 15°S and 15°N) of (a) the specific humidity anomalies between ERA5 and ERA-Interim reanalyses, (b) E6_LR_ERA5 and E6_LR_ERAI simulations, and of (c) the differences in modeled δ 18 O v between E6_LR_ERA5 and E6_LR_ERAI simulations.
are found over the South American continent (70°W) at 900 hPa level. Slightly positive anomalies are found over the Atlantic Ocean at 600 hPa level (less than 2 × 10 −4 kg kg −1 ) but not as strong as in the reanalyses anomalies (between 2 and 5 × 10 −4 kg kg −1 ), and only modeled negative anomalies are found over Africa. These changes in the distribution of modeled water vapor and of its transport have impacts on the modeled δ 18 O v values, mostly at pressure levels from 500 to 200 hPa (Figure 9c). Slightly negative anomalies (less than −0.50‰) are found over the Pacific, Atlantic and Indo-Pacific areas. Stronger positive δ 18 O v anomalies are located over the eastern coast of Brazil at 55°W in the high troposphere (∼1.0‰) and spread down to 750 hPa, probably related to the stronger transport of marine water vapor (i.e., enriched in heavy isotopes) to Brazil due to nudging to ERA5 instead of ERA-Interim. As a counterpart, negative δ 18 O v anomalies are located at 80°W in the western part of the South American continent, consistent with the increase in precipitation from E6_LR_ERAI to E6_LR_ERA5 (Figure 7c) depleting the remaining water vapor in heavy isotopes. Further investigations by comparing these findings with satellite-based observations of stable water isotopes could be performed in the future. A similar positive/negative anomalies dipole pattern in δ 18 O is modeled in precipitation (Figure 7a), confirming the influence of water vapor transport change on δ 18 O p between the E6_LR_ERAI and E6_LR_ERA5 simulations.

Conclusions
In this study we have presented new isotope-enabled ECHAM6-wiso simulations nudged to ERA5 reanalyses for the period 1979-2018 and investigated the effects of the changes in the reanalyses data used for the nudging, general improvements of the ECHAM6 model as compared to its predecessor and the spatial model resolution on the isotope-related variables. Since the model version published in Cauquoin et al. (2019), ECHAM6-wiso has been updated in several aspects to be consistent with the more recent isotope observation campaigns (Bonne et al., 2019). The kinetic fractionation factors for the evaporation over the ocean are now assumed as independent of wind speed. A similar model-data agreement is obtained compared to the model results using the standard model approach by Merlivat and Jouzel (1979). Moreover, the isotopic content of snow on sea ice is now taken into for sublimation processes in sea ice covered regions. This leads to a stronger depletion of surface water vapor over such sea ice covered areas. Finally, the supersaturation equation has been slightly re-tuned to better reproduce the isotope delta values measured in Antarctic snow. All the simulations show a reasonably good agreement with the various selected observational groundbased isotope data sets. On the global scale, our modeled water isotope δ values in precipitation and surface water vapor are rather similar for the different ERA reanalyses data used for the nudging, even if significant changes in tropical rainfall patterns and lower surface temperatures over land can be observed.
Despite small apparent variations in surface variables by switching from ERA-Interim to ERA5 reanalyses for the nudging, changes in tropical water vapor distribution, like increased transport to the Amazonian area, lead to enriched water vapor in heavy isotopes in the high troposphere over the eastern coast of Brazil and depleted water vapor in heavy isotopes over the western coast of South America. The modeled specific humidity anomalies between the two ECHAM6-wiso simulations nudged to ERA-Interim and ERA5 have some similarities with the anomalies from the ERA reanalyses, like increased specific humidity over Amazonia, negative anomalies over Southern Sahara in the low troposphere and uniform negative anomalies at ∼800 hPa pressure level. On the other hand, positive anomalies in specific humidity between ERA reanalyses over the Atlantic Ocean, Africa and Indo-Pacific regions higher in the troposphere are not reflected in the ECHAM6-wiso anomalies. The amplitude of the changes in modeled specific humidity between the ECHAM6-wiso simulations nudged to ERA-Interim and ERA5 is also smaller than the one found in the ERA reanalyses anomalies. These differences are most likely explained by different formulations of atmospheric hydrological processes (e.g., evaporation, condensation, cloud formation) between the ERA and ECHAM6-wiso. The differences in specific humidity are then imprinted in the spatial distribution of water vapor isotopes in the atmosphere, too.
It is important to note that our findings about the switching from ERA-Interim to ERA5 reanalyses for the nudging are not valid with ECHAM6-wiso, only. Very comparable results were obtained with ECHAM5wiso, too (see Supporting Information S1 for details). As the two ECHAM versions are quite different in many physical model aspects related to the radiation scheme, the cloud cover, or the land surface processes (Stevens et al., 2013), we rate our results as robust. Comparable studies with other isotope-enabled general circulation models are necessary to quantify the general validity of our results more firmly.
The changes of modeled water vapor transport over the tropics when using ERA5 for nudging instead of ERA-Interim is an important aspect for studies using isotopic records in this region, like tree rings or speleothem archives in South America, or corals in tropical latitudes. The improved representation of extreme events like tropical cyclones in ERA5 (Hersbach et al., 2020) is also a great opportunity to study the key processes controlling the isotopic variability in precipitation and water vapor during such events, with potential applications for the interpretation of paleo-hydroclimate across the tropics (Sánchez-Murillo et al., 2019).
Some substantial changes are found in modeled isotope values at the surface between the predecessor ECHAM5-wiso and the current model version ECHAM6-wiso, even if both models are nudged to the same ERA5 reanalyses data. These changes can be related to various general ECHAM model updates of key physical processes, for example, the radiation scheme, land-surface processes, or the chosen vertical discretization (Stevens et al., 2013). It is worthwhile to note that the general ECHAM6 model improvements do not always automatically lead to an improved simulation of related isotope signals.
With respect to the spatial model resolution, improvements are seen for the higher spatial resolution simulation ECHAM6-wiso T127L95 for annual mean values and temporal variations of isotopic δ values at seasonal and (sub-)daily time scales. Thus, we conclude that especially this new ECHAM6-wiso simulation with high spatial resolution and nudged to the most recent ERA5 data set will be extremely valuable for follow-up studies involving water isotopes measured in various natural archives as well as for the analysis of water vapor transport variability in the tropics or/ and in the higher atmosphere. The extension of the ERA5 data set back to the year 1950 in the near future and the continuous production of new reanalysis data with a delay of 2-3 months will make it possible to have high resolution ECHAM6-wiso simulation products on the entire period from 1950 to present time, soon.

Data Availability Statement
The ERA5 reanalyses for the nudging have been downloaded from the German Climate Computing Center (DKRZ), generated using Copernicus Climate Change Service Information (2020). The ECHAM5-wiso and ECHAM6-wiso simulations have been performed at the Alfred Wegener Institute (AWI) supercomputing centre. The ECHAM model code is available under a version of the MPI-M software license agreement (https://www.mpimet.mpg.de/en/science/models/license/). The code of the isotopic version ECHAM6wiso is available upon request on the AWI's GitLab repository (https://gitlab.awi.de/mwerner/mpi-esmwiso, Cauquoin et al., 2019). Observational datasets for this research are included in these papers ( 2008) and Wei et al. (2015Wei et al. ( , 2016. GNIP data are available at https://nucleus.iaea. org/wiser. The model data used in this study are available on the Zenodo database: https://doi.org/10.5281/ zenodo.5636328.