From climatological to small-scale applications: simulating water isotopologues with ICON-ART-Iso (version 2.3)

We present the new isotope-enabled model ICON-ART-Iso. The physics package of the global ICOsahedral Nonhydrostatic (ICON) modeling framework has been extended to simulate passive moisture tracers and the stable isotopologues HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$. The extension builds on the infrastructure provided by ICON-ART, which allows for high flexibility with respect to the number of related water tracers that are simulated. The physics of isotopologue fractionation follow the model COSMOiso. We first present a detailed description of the physics of fractionation that have been implemented in the model. The model is then evaluated on a range of temporal scales by comparing with measurements of precipitation and vapor.

A multi-annual simulation is compared to observations of the isotopologues in precipitation taken from the station network GNIP (Global Network for Isotopes in Precipitation). ICON-ART-Iso is able to simulate the main features of the seasonal cycles in δD and δ¹⁸O as observed at the GNIP stations. In a comparison with IASI satellite retrievals, the seasonal and daily cycles in the isotopologue content of vapor are examined for different regions in the free troposphere. On a small spatial and temporal scale, ICON-ART-Iso is used to simulate the period of two flights of the IAGOS-CARIBIC aircraft in September 2010, which sampled air in the tropopause region influenced by Hurricane Igor. The general features of this sample as well as those of all tropical data available from IAGOS-CARIBIC are captured by the model.

The study demonstrates that ICON-ART-Iso is a flexible tool to analyze the water cycle of ICON. It is capable of simulating tagged water as well as the isotopologues HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$.

1 Introduction

Water in gas, liquid and frozen form is an important component of the climate system. Ice caps and snow-covered surfaces strongly influence the albedo of the surface (Kraus, 2004), the oceans are unmatched water reservoirs, which dissolve trace substances (Jacob, 1999) and redistribute heat (Pinet, 1993), and all animal and plant life depends on liquid water. The atmosphere is by mass the smallest compartment of the hydrological cycle, but it is this compartment that serves to transfer water between the spheres of liquid, frozen and biologically bound water on the earth's surface (Gat, 1996). For atmospheric processes themselves, water is also of great importance. It is the strongest greenhouse gas (Schmidt et al., 2010) and distributes energy through the release of latent heat (Holton and Hakim, 2013), while liquid and frozen particles influence the radiative balance (Shine and Sinha, 1991).

A correct description of the atmospheric water cycle is therefore necessary for the understanding and simulation of the atmosphere and the climate system (Riese et al., 2012; Sherwood et al., 2014). The stable isotopologues of water are unique diagnostic tracers that provide a deeper insight into the water cycle (Galewsky et al., 2016). Because of the larger molar mass of the heavy isotopologues, their ratio to (standard) water is changed by phase transitions. This change in the ratio is termed fractionation. Considering the isotopologue ratio of the heavy isotopologues in vapor and precipitation (liquid or ice) provides an opportunity to develop an advanced understanding of the processes that shape the water cycle.

Pioneering research on measuring the heavy isotopologues of water started in the 1950s and first examined the isotopologues in precipitation (Dansgaard, 1954, 1964). Theoretical advances on the microphysics (Jouzel and Merlivat, 1984; Jouzel et al., 1975) and surface evaporation (Craig and Gordon, 1965) enabled the implementation of heavy isotopologues in global climate models (Joussaume and Jouzel, 1993; Joussaume et al., 1984). Since then, measurement techniques and modeling of the isotopologues have advanced. Measurements of the isotopic content of vapor first required cryogenic samplers (Dansgaard, 1954), but in the last 15 years laser absorption spectroscopy has made in situ observations possible (Dyroff et al., 2010; Lee et al., 2005). Today, the isotopologue content in atmospheric vapor can also be derived from satellite measurements (Gunson et al., 1996; Schneider and Hase, 2011; Steinwagner et al., 2007; Worden et al., 2006). Many global and regional circulation models have been equipped to simulate the atmospheric isotopologue distribution, focusing on global-scale (Risi et al., 2010; Werner et al., 2011) or regional phenomena (Blossey et al., 2010; Pfahl et al., 2012). Despite this progress, the potential of isotopologues to improve the understanding and physical description of single processes “remains largely unexplored” (Galewsky et al., 2016). A more extensive literature overview on the subject is given by Galewsky et al. (2016).

We present ICON-ART-Iso, the newly developed, isotopologue-enabled version of the global ICOsahedral Nonhydrostatic (ICON) modeling framework (Zängl et al., 2015). By design, ICON is a flexible model capable of simulations from climatological down to turbulent scales (Heinze et al., 2017). The advection scheme of ICON has been designed to be mass conserving (Zängl et al., 2015), which is essential for the simulation of water isotopologues (Risi et al., 2010). ICON-ART-Iso builds on the flexible infrastructure provided by the extension ICON-ART (Rieger et al., 2015; Schröter et al., 2018), which has been developed to simulate aerosols and trace gases.

By equipping ICON with the capabilities to simulate water isotopologues, a first step is made to a deeper understanding of the water cycle. From the multitude of isotopologue-enabled global models (see Galewsky et al., 2016 for an overview), ICON-ART-Iso stands out because of its non-hydrostatic base model core, enabling simulations with fine horizontal resolution on a global grid. Its flexible design allows for the simulation of diagnostic evaporation tracers as well as the isotopologues HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ during a single simulation.

This article first gives some technical details on ICON and ICON-ART. This is followed by a detailed description of the physics special to ICON-ART-Iso, which have been implemented in ICON to simulate isotopologues (Sect. 2).

The remaining sections describe model results and first validation studies: Sect. 3.1 looks at passive moisture tracers. Focus is laid on the source regions – ocean or land – of the water that forms precipitation. The next section (Sect. 3.2) compares data from a simulation spanning more than 10 years on a coarse grid to measurements from different stations of the GNIP network. A further validation with measurements is performed in Sect. 3.3. Retrievals from IASI satellite measurements are compared with ICON-ART-Iso results for 2 weeks in winter and summer 2014, considering the seasonal and daily cycle in different regions. Section 3.4 then discusses the comparison with IAGOS-CARIBIC measurements. In situ data from two flights are compared with results of ICON-ART-Iso simulations. Section 4 summarizes and concludes the study.

2 The model ICON-ART-Iso

This section presents the technical and physical background of the model ICON-ART-Iso. First, ICON and the extension ICON-ART are introduced. Next, general thoughts on simulating a diagnostic water cycle are presented. Starting in Sect. 2.3, the main processes that influence the distribution of the isotopologues are discussed in separate sections: surface evaporation, saturation adjustment, cloud microphysics and convection. To close this technical part, Sect. 2.7 discusses the initialization of the model.

2.1 Introduction to the modeling framework ICON-ART

ICON-ART-Iso is the isotope-enabled version of the model ICON. ICON is a new non-hydrostatic general circulation model developed and maintained in a joint effort by the Deutscher Wetterdienst (DWD) and Max Planck Institute for Meteorology (MPI-M). Its horizontally unstructured grid can be refined locally by one-way or two-way nested domains with a higher resolution. The model is applicable from global to turbulent scales: at DWD, ICON is used operationally for global numerical weather prediction (currently 13 km horizontal resolution, with a nest of 6.5 km resolution over Europe). Klocke et al. (2017) show the potential of using ICON for convection-permitting simulations and it already proved successful as a Large eddy simulation (LES) model (Heinze et al., 2017). It is currently also being prepared for climate projections at MPI-M. More details on ICON are given by Zängl et al. (2015).

ICON-ART-Iso builds on the numerical weather prediction physics parameterization package of ICON. The physical parameterizations that have been implemented for the simulation of isotopologues mainly correspond to those of the model COSMOiso as presented by Pfahl et al. (2012). As the same parameterizations have been described before, the following subsections give only a short summary of each of the different fractionation processes.

In ICON, all tracer constituents are given as mass fractions $q_x=\frac{{\mathit{\rho }}_x}{\mathit{\rho }}$, where $\mathit{\rho }={\sum }_x{\mathit{\rho }}_x$ is the total density, including all water constituents x. To discriminate the values of heavy isotopologues, these will be denoted by the index h, while standard (light) water will be indexed by the letter l. ICON standard water is identified with the light isotopologue ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$, which is a very good assumption also made by Blossey et al. (2010) and Pfahl et al. (2012): standard water is much more abundant than the lighter isotopologues, with a ratio of 1 to $\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{4}}$ for HDO and $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{3}}$ for ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ (Gonfiantini et al., 1993). Water in ICON-ART-Iso exists in seven different forms (vapor, cloud water, ice, rain, snow, graupel and hail), each of which is represented by one tracer for standard water and an additional tracer for each of the isotopologues. The amount of the isotopologues is expressed relative to standard water by the isotopologue ratio $R={}^{\mathrm{h}}{q}_x/{}^{\mathrm{l}}{q}_x$. This is referenced to standard ratios of the Vienna Mean Ocean Standard Water (R_VSMOW) in the δ notation: $\mathit{\delta }=R_{\text{sample}}/R_{\text{VSMOW}}-\mathrm{1}$, with δ values then given in per mil. If not noted otherwise, δ values are always evaluated for the vapor phase in this paper, which is why this specification is omitted throughout the text.

As in the current version of ICON-ART (Schröter et al., 2018), an XML table is used to define the settings for each of the isotopologues. While this paper mostly discusses realizations of HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$, this choice is technically arbitrary. The XML table is used to define the tracers at runtime, making a recompilation of the model unnecessary. All tuning parameters can be specified separately for each isotopologue in the XML table and the number of realizations is limited only by the computational resources. Each parameterization describing fractionation can also be turned off separately for each isotopologue, making very different experiments possible during one simulation. This makes the model very flexible and allows for the use of several different water tracers during one model run.

2.2 Simulating a diagnostic water cycle

The isotopologues are affected by all the processes that also influence standard water in ICON: surface evaporation, saturation adjustment to form clouds, cloud microphysics and convection. Each of these main processes is represented by several parameterizations. Some of these parameterizations include phase changes of or to vapor and in turn lead to a change in the isotopologue ratio, which is termed isotopic fractionation. In addition, advection and turbulent diffusion are non-fractionating processes that change the spatial distribution of all trace substances.

An important prerequisite to a simulation of water isotopologues is a good implementation of advection. ICON-ART makes use of the same numerical methods that are used for advecting the hydrometeors in ICON itself. These ensure local mass conservation (Zängl et al., 2015) and mass-consistent transport. The latter is achieved by making use of the same mass flux in the discretized continuity equations for total density and partial densities, respectively (Lauritzen et al., 2014). The advection schemes implemented in ICON conserve linear correlations between tracers and ensure the monotonicity of each advected tracer. Note, however, that this does not guarantee monotonicity of the isotopologue ratios (Morrison et al., 2016).

The parameterizations that influence the water cycle also include processes that do not fractionate. For all non-fractionating processes, the transfer rate ^hS of the heavier isotopologues is defined by Eq. (1).

$$\begin{array}{}\text{(1)}& {}^{\mathrm{h}}S={}^{\mathrm{l}}S\cdot R_{\text{source}}\end{array}$$

Here, ^lS is the transfer rate of ICON standard water, while R_source is the isotopologue ratio in the source reservoir of the transfer.

In order to turn any fractionating processes into a non-fractionating one, its respective equation for the transfer rate of the heavy isotopologues can be replaced with Eq. (1). This has been implemented as an option in all processes that describe fractionation, which are explained below. If all processes are set to be non-fractionating in this way, the isotopologue ratio does not change and the species will resemble the standard water of ICON. This is an important feature that can be used to test the model for self-consistency or to investigate source regions with diagnostic moisture tracers, so-called tagged water (Bosilovich and Schubert, 2002). An application of this will be shown in Sect. 3.1.

Whenever phase changes occur that include the vapor phase, the isotopologue ratio changes because the heavier isotopologues have different diffusion constants and a different saturation vapor pressure compared to standard water. For the diffusion constant ratio, two choices have been implemented for HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$, making available the values of Merlivat and Jouzel (1979) or Cappa et al. (2003). The differences in saturation pressure are expressed by the equilibrium fractionation factor α, which is the ratio of isotopologue ratios in thermodynamic equilibrium (Mook, 2001); see Eq. (2).

$$\begin{array}{}\text{(2)}& \mathit{\alpha }={\displaystyle \frac{R_{\mathrm{v}}}{R_{\text{cond}}}}<\mathrm{1}\end{array}$$

Here, R_v stands for the isotopologue ratio in the vapor phase, while R_cond stands for that in the condensed phase. The ratio α depends on temperature and is different over water and over ice (termed α_liq and α_ice). The parameterizations by Majoube (1971) and by Horita and Wesolowski (1994) have been implemented for α_liq and those by Merlivat and Nief (1967) for α_ice. Note the definition for α given in Eq. (2) is also used in COSMOiso (Pfahl et al., 2012) and is the inverse of the definition used by others, e.g., by Blossey et al. (2010).

2.3 Surface evaporation

Surface evaporation is the source for the atmospheric water cycle. In ICON-ART-Iso, the evaporative surface flux is split into evaporation from land and water surfaces, transpiration from plants, and dew and rime formation. Transpiration is considered a non-fractionating process (Eq. 1), which is an assumption also made by Werner et al. (2011) and Pfahl et al. (2012). Dew and rime formation (and condensation on the ocean surface) are considered to fractionate according to equilibrium fractionation (Eq. 2). For the evaporation part of the full surface flux, two parameterizations have been implemented (Merlivat and Jouzel, 1979; Pfahl and Wernli, 2009). Both build on the Craig–Gordon model (Craig and Gordon, 1965; Gat, 2010). Equation (3) gives the general expression for R_evap.

$$\begin{array}{}\text{(3)}& R_{\text{evap}}=k\cdot {\displaystyle \frac{{\mathit{\alpha }}_{\text{liq}}{R}_{\text{surf}}-h{R}_{\mathrm{v}}}{\mathrm{1}-h}}\end{array}$$

Here, h is the specific humidity of the lowest model layer relative to the specific humidity at the surface and k is the nonequilibrium fractionation factor. The two parameterizations differ in their description of k. While Merlivat and Jouzel (1979) give a parameterization that depends on the surface wind, Pfahl and Wernli (2009) have simplified this to be wind speed independent. In summary, Eq. (4) is used to calculate the surface flux of the isotopologues, ^hF^tot.

$$\begin{array}{ll}{}^{\mathrm{h}}{F}^{\mathrm{tot}}=& {\displaystyle }{}^{\mathrm{l}}{F}^{\text{evap}}\cdot R_{\text{evap}}+{}^{\mathrm{l}}{F}^{\text{transp}}\cdot R_{\text{surf}}+{}^{\mathrm{l}}{F}^{\text{dew}}\\ \text{(4)}& {\displaystyle }& {\displaystyle }\cdot {\displaystyle \frac{R_{\mathrm{v}}}{{\mathit{\alpha }}_{\text{liq}}}}+{}^{\mathrm{l}}{F}^{\text{rime}}\cdot {\displaystyle \frac{R_{\mathrm{v}}}{{\mathit{\alpha }}_{\text{ice}}}}\end{array}$$

For transpiration and evaporation, the isotopologue ratio of the surface water and groundwater (R_surf) is necessary. The surface model TERRA (included in ICON) was not extended with isotopologues, so R_surf is not available as a prognostic variable. Over land, it is therefore approximated by R_VSMOW in Eqs. (3) and (4). Of course, this is a simplification that allows for the testing of the atmospheric physics package and will be developed further. Over the ocean, the dataset provided by LeGrande and Schmidt (2006) has been implemented. Values for HDO are given in this dataset, while those for ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ are determined from the relationship given by the global meteoritic water line (GMWL), $\mathit{\delta }\text{D}=\mathrm{8}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{10}$ ‰ (Craig, 1961).

2.4 Saturation adjustment

Cloud water is formed by saturation adjustment in ICON. Vapor in excess of saturation vapor pressure is transferred to cloud water and temperature is adjusted accordingly. This is repeated in an iterative procedure. For the isotopologues, the iteration does not have to be repeated. Instead, Eq. (5) is applied directly using the adjusted values of ICON water. This is the same equation used in COSMOiso (Pfahl et al., 2012) and by Blossey et al. (2010).

$$\begin{array}{}\text{(5)}& {}^{\mathrm{h}}{q}_{\mathrm{c}}={\displaystyle \frac{{}^{\mathrm{h}}{q}_{\mathrm{v}}+{}^{\mathrm{h}}{q}_{\mathrm{c}}}{\mathrm{1}+{\mathit{\alpha }}_{\text{liq}}\frac{{}^{\mathrm{l}}{q}_{\mathrm{v}}}{{}^{\mathrm{l}}{q}_{\mathrm{c}}}}}\end{array}$$

2.5 Microphysics

Several grid-scale microphysical schemes are available in ICON. ICON-ART-Iso makes use of the two-moment scheme by Seifert and Beheng (2006). This scheme computes the mass and number densities of vapor, cloud water, rain and four ice classes (ice, snow, graupel and hail) and can be used to simulate aerosol–cloud interaction; see Rieger et al. (2017). As the isotopologues are diagnostic values, the number densities do not have to be simulated separately. The two-moment scheme describes more than 60 different processes, but only those processes that include the vapor phase lead to fractionation. All others are described by Eq. (1) in the model. Isotopic effects also occur during freezing of the liquid phase (Souchez and Jouzel, 1984; Souchez et al., 2000), but this is neglected due to the low diffusivities, as in COSMOiso. In accordance with Blossey et al. (2010) and Pfahl et al. (2012), sublimation is also assumed not to fractionate. Condensation to form liquid water happens only during the formation of cloud water and is accounted for by the saturation adjustment. The fractionating processes that remain are ice formation by nucleation, vapor deposition (on all four ice classes) and evaporation of liquid hydrometeors. Besides rain, a fraction of the three larger ice classes (snow, graupel, hail) can evaporate after melting. This liquid water fraction is currently not a prognostic variable.

The two-moment scheme by Seifert and Beheng (2006) uses mass densities instead of mass ratios, so we adopt the change in notation here, denoting mass densities by ρ. Vapor pressures are denoted by e. The star (^*) indicates values at saturation with respect to liquid (index l) or ice (index i).

For evaporation of rain and melting hydrometeors, the semiempirical parameterization of Stewart (1975) has been implemented and is discussed in this paper. It allows for the exchange of heavy isotopologues with the surroundings in supersaturated as well as subsaturated conditions. The corresponding transfer rate is given in Eq. (6). The equation is given in the formulation for the evaporation of rain, with details on the evaporation of melting ice hydrometeors explained below. In this process, it is assumed that the isotopic content within each droplet is well mixed, which is a simplification when compared to, e.g., Lee and Fung (2008).

$$\begin{array}{}\text{(6)}& {\displaystyle }& {}^{\mathrm{h}}{S}_{\mathrm{r}}^{\text{evap}}=A{\left({\displaystyle \frac{{}^{\mathrm{h}}D}{{}^{\mathrm{l}}D}}\right)}^n\left[R_{\mathrm{r}}{\mathit{\alpha }}_{\text{liq}}^{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l},\mathrm{\infty }}^{*}-{}^{\mathrm{h}}{\mathit{\rho }}_{\mathrm{v}}\right]\text{(7)}& {\displaystyle }& {\displaystyle }A={\displaystyle \frac{\mathrm{4}\mathit{\pi }a{}^{\mathrm{l}}f{}^{\mathrm{l}}D}{\mathrm{1}+B_{\mathrm{l}}}}\text{(8)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{l}}={\displaystyle \frac{{}^{\mathrm{l}}D{L}_{\mathrm{e}}^{\mathrm{2}}{}^{\mathrm{l}}{e}_{\mathrm{l},\mathrm{\infty }}^{*}}{k_{\mathrm{a}}{\mathcal{R}}_{\mathrm{v}}^{\mathrm{2}}{T}_{\mathrm{\infty }}^{\mathrm{3}}}}\end{array}$$

Here, a is the radius of the hydrometeor, ^lf is the ventilation factor, R_v and R_r are the isotopologue ratios in the vapor and the hydrometeor, ℛ_v is the gas constant of water vapor, and L_e and k_a are the latent heat of evaporation and the heat conductivity in air. The index ∞ indicates that values are evaluated for the surroundings. The ratio of the diffusion constants D is given by the literature values cited above and can be chosen at the time of simulation for each isotopologue. The tuning parameter n is set to 0.58 by default (Stewart, 1975), but can be changed at runtime.

Note that an alternative parameterization to describe the fractionation of evaporating or equilibrating hydrometeors (that of Blossey et al., 2010) has also been implemented in the model. For completeness, the physics of this parameterization are briefly explained in Appendix B in comparison to Stewart (1975). An investigation of this parameterization and the difference to Stewart (1975) will be provided in a later study.

The underlying equation for both parameterizations is derived from the fundamentals of cloud microphysics (Pruppacher and Klett, 2012). It is also used in the microphysical scheme of ICON, in which ${}^{\mathrm{l}}{S}_x^{\text{evap}}=A({}^{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l},\mathrm{\infty }}^{*}-{}^{\mathrm{l}}{\mathit{\rho }}_{\mathrm{v}})$. The definitions of a and ^lf_v depend on whether S_x is calculated for rain (Seifert, 2008) or melting ice class hydrometeors (Seifert and Beheng, 2006). For melting ice class hydrometeors, the melting temperature of ice (T₀=273.15 K) is used for the calculation of α_liq and in place of T_∞. This implies an additional factor of T_∞∕T₀ for melting ice hydrometeors, as ${}^{\mathrm{l}}{e}_{\mathrm{l},\mathrm{\infty }}^{*}$ is always evaluated at T_∞. Equation (6) otherwise also holds true for the evaporation of melting ice hydrometeors.

Fractionation during the nucleation of ice particles or deposition on one of the four ice class hydrometeors is parameterized following Blossey et al. (2010), as in COSMOiso (Pfahl et al., 2012). The flux is assumed to interact only with the outermost layer of the hydrometeor, the isotopologue ratio of which is set to be identical to that of the depositional flux. The transfer rate ${}^{\mathrm{h}}{S}_x^{\text{ice}}$ is then given by Eq. (9) with the fractionation factor α_k as given in Eq. (10). All symbols are used as above, with L_s being the latent heat of sublimation.

$$\begin{array}{}\text{(9)}& {\displaystyle }& {}^{\mathrm{h}}{S}_x^{\text{ice}}={\mathit{\alpha }}_k{R}_{\mathrm{v}}^{\mathrm{l}}{S}_x^{\text{ice}}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_k={\displaystyle \frac{\left(\mathrm{1}+B_{\mathrm{i}}\right){}^{\mathrm{l}}{\mathcal{S}}_{\mathrm{i}}}{\frac{{}^{\mathrm{l}}f{}^{\mathrm{l}}D}{{}^{\mathrm{h}}f{}^{\mathrm{h}}D}\left({}^{\mathrm{l}}{\mathcal{S}}_{\mathrm{i}}-\mathrm{1}\right)+{\mathit{\alpha }}_{\text{ice}}\left(\mathrm{1}+B_{\mathrm{i}}{}^{\mathrm{l}}{\mathcal{S}}_{\mathrm{i}}\right)}}\text{(11)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{i}}={\displaystyle \frac{{}^{\mathrm{l}}{D}_{\mathrm{v}}{L}_s^{\mathrm{2}}{}^{\mathrm{l}}{e}_{i,\mathrm{\infty }}^{*}}{k_{\mathrm{a}}{\mathcal{R}}_{\mathrm{v}}^{\mathrm{2}}{T}_{\mathrm{\infty }}^{\mathrm{3}}}}\end{array}$$

2.6 Convection

ICON uses the Tiedtke–Bechtold scheme for simulating convective processes (Bechtold et al., 2014; Tiedtke, 1989). The scheme uses a simple cloud model considering a liquid fraction in cloud water (denoted here by ω) and the remaining solid fraction (1ω). Fractionation happens during convective saturation adjustment (during initialization of convection and in updrafts), in saturated downdrafts and in evaporation below cloud base. The parameterizations are the same as those implemented by Pfahl et al. (2012) in COSMOiso.

Convective saturation adjustment calculates equilibration between vapor and the total condensed water (liquid and ice). The parameterization used for grid-scale adjustment therefore has to be expanded in order to be used in convection if the liquid water fraction is smaller than one. The isotopologue ratio is determined over liquid and ice particles separately. A closed system approach (Gat, 1996) is used for the liquid fraction ($R_{\mathrm{v}}^{\text{by liq}}$ of Eq. 12). The underlying assumption for Eq. (13) used for the ice fraction is a Rayleigh process with the kinetic fractionation factor α_eff following Jouzel and Merlivat (1984). The two are then recombined according to the fraction of liquid water, following Eq. (14). This procedure has been adopted from COSMOiso (Pfahl et al., 2012).

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{v}}^{\text{by liq}}=R_{\mathrm{v}}^{\text{old}}{\displaystyle \frac{{\mathit{\alpha }}_{\text{liq}}}{\mathrm{1}+\frac{{}^{\mathrm{l}}{q}_{\mathrm{v}}^{\text{new}}}{{}^{\mathrm{l}}{q}_{\mathrm{v}}^{\text{old}}}({\mathit{\alpha }}_{\text{liq}}-\mathrm{1})}}\text{(13)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{v}}^{\text{by ice}}=R_{\mathrm{v}}^{\text{old}}{\left({\displaystyle \frac{{}^{\mathrm{l}}{q}_{\mathrm{v}}^{\text{new}}}{{}^{\mathrm{l}}{q}_{\mathrm{v}}^{\text{old}}}}\right)}^{{\mathit{\alpha }}_{\text{eff}}-\mathrm{1}}\text{(14)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{v}}=(\mathrm{1}-\mathit{\omega })\cdot R_{\mathrm{v}}^{\text{by ice}}+\mathit{\omega }\cdot R_{\mathrm{v}}^{\text{by liq}}\end{array}$$

Here, the indices “old” and “new” denote the values of the respective variables before and after the convective saturation adjustment. The factor α_eff, which appears in Eq. (13), is determined by Eq. (15). The supersaturation with respect to ice, ξ_ice, is calculated from Eq. (16), where T₀=273.15 K is used. The tuning parameter λ is set to 0.004 in the standard setup, following Pfahl et al. (2012) and Risi et al. (2010).

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\text{eff}}={\displaystyle \frac{{\mathit{\xi }}_{\text{ice}}\mathit{\zeta }}{{\mathit{\xi }}_{\text{ice}}-\mathrm{1}+{\mathit{\alpha }}_{\text{ice}}\mathit{\zeta }}}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\xi }}_{\text{ice}}=\mathrm{1}-\mathit{\lambda }(T-T_{\mathrm{0}})\end{array}$$

Convective downdrafts are assumed to remain saturated by continuously evaporating precipitation (Tiedtke, 1989). In these saturated downdrafts, equilibrium fractionation is applied for the liquid fraction, while the ice fraction is assumed to sublimate without fractionation.

Evaporation of precipitation below cloud base is an important process for several reasons: it leads to a drop in the temperature and therefore influences dynamics, but is also important for the isotopic composition (Risi et al., 2008). To describe fractionation here, the parameterization by Stewart (1975) is again applied to the liquid fraction. Different to Eq. (6) for evaporation during microphysics, the integrated form is now applied. Following Stewart (1975), the ratio in the liquid part of the general hydrometeor after evaporation $R_{\text{adj}}^{\text{liq}}$ is given with Eq. (17). Here, f is the fraction of remaining condensate. $R_{\text{hyd}}^{\text{old}}$ is the isotopologue ratio in the hydrometeor before adjustment and RH is the relative humidity calculated as the vapor pressure over saturation vapor pressure.

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }{R}_{\text{adj}}^{\text{liq}}=\mathit{\gamma }{R}_{\mathrm{v}}+f^{\mathit{\beta }}\left(R_{\text{hyd}}^{\text{old}}-\mathit{\gamma }{R}_{\mathrm{v}}\right)\text{(18)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{\text{RH}}{{\mathit{\alpha }}_{\text{liq}}-\mathit{\mu }}}\text{(19)}& {\displaystyle }& {\displaystyle }\mathit{\beta }={\displaystyle \frac{{\mathit{\alpha }}_{\text{liq}}-\mathit{\mu }}{\mathit{\mu }}}\text{(20)}& {\displaystyle }& {\displaystyle }\mathit{\mu }=\left(\mathrm{1}-\text{RH}\right){\left({\displaystyle \frac{{}^{\mathrm{h}}D}{{}^{\mathrm{l}}D}}\right)}^{-n}\end{array}$$

Using Eq. (17), the isotopologue ratio in the adjusted hydrometeor is given with Eq. (21). The ice fraction is assumed to sublimate without fractionation, maintaining its isotopologue ratio.

$$\begin{array}{}\text{(21)}& R_{\text{adj}}=(\mathrm{1}-\mathit{\omega })R_{\text{hyd}}^{\text{old}}+\mathit{\omega }{R}_{\text{adj}}^{\text{liq}}\end{array}$$

Following Pfahl et al. (2012), an additional equilibration has been implemented to determine the final isotopologue ratio of the hydrometeors, which is given in Eq. (22). The parameter ξ_add is a tuning parameter that is set to 0.5 in the standard setup.

$$\begin{array}{}\text{(22)}& R_{\text{adj}}^{\text{final}}=R_{\text{adj}}+{\mathit{\xi }}_{\text{add}}\cdot \mathit{\omega }\left({\displaystyle \frac{R_{\mathrm{v}}}{{\mathit{\alpha }}_{\text{liq}}}}-R_{\text{adj}}\right)\end{array}$$

2.7 Initialization of the isotopologues

A meaningful initialization is an important prerequisite for any simulation, also of the isotopologues. In addition to an initialization with a constant ratio to standard water, the isotopologues can be initialized with the help of mean measured δ values. Values at the lowest model level, the tropopause level (Holton et al., 1995) and the model top are prescribed for vapor, and linear and log-linear interpolation is applied below and above the tropopause, respectively. Values for the tropopause level and the model top are taken from MIPAS measurements (Steinwagner et al., 2007), and the value at the lowest level is a standard value taken from Gat (2010). All values are given in Table 1. By using the local tropopause height, an adaptation to the local meteorological situation is ensured. To calculate the δ value of the hydrometeors, a constant offset is applied to the local δ value of vapor. The literature provides values for HDO, while those for ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ are determined from the relationship given by the global meteoritic water line (Craig, 1961).

Gat (2010)Steinwagner et al. (2007)Steinwagner et al. (2007)Gat (2010)

Table 1 Values for the initialization with mean measured δ values. The literature provides the values for HDO, and values for ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ have been determined from GMWL (Craig, 1961).

3 Model evaluation results

In the following sections, we present the first results and comparisons of model simulations with measurements spanning several spatiotemporal scales: Sect. 3.1 shows how the model-simulated diagnostic H₂O can be used to investigate source regions of the (modeled) water cycle. Section 3.2 compares results for precipitation from the same long model integration with measurements taken from the GNIP network (IAEA/WMO, 2017; Terzer et al., 2013). Section 3.3 looks at seasonal and regional differences by comparing model output with pairs of {H₂O,δD} derived from IASI satellite measurements (Schneider et al., 2016). Finally, Sect. 3.4 presents a first case study, in which simulated values of δD are compared with measurements from the IAGOS-CARIBIC project (Brenninkmeijer et al., 2007). All simulations discussed here are free-running.

In the following sections, we focus on H₂O (ICON standard water), HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$. The settings for each isotopologue are defined at runtime, which is why the specifications for the simulations are given here. The diffusion constant ratio is set to the values of Merlivat and Jouzel (1979) and the equilibrium fractionation is parameterized following Majoube (1971) over liquid water and following Merlivat and Nief (1967) over ice. Surface evaporation is described by the parameterization of Pfahl and Wernli (2009). The parameterization by Merlivat and Jouzel (1979) has little influence on the values in the free troposphere and is not discussed further. The dataset by LeGrande and Schmidt (2006) for the isotopic content of the ocean surface is used for all isotopologues. The grid-scale evaporation of hydrometeors is described by the parameterization by Stewart (1975) if not noted differently.

In addition to ICON standard water, three diagnostic sets of water tracers are simulated. All fractionation is turned off, so they resemble H₂O. But the evaporation and initialization are different: water indexed by init (as in q^init) is set to ICON water at initialization, but evaporation is turned off. In the course of the simulation, water of this type precipitates out of the model atmosphere. Water indexed by ocn and lnd is initialized with zero and evaporates from the ocean (q^ocn) and land areas (q^lnd), respectively. The sum of q^init, q^ocn and q^lnd always equals the mass mixing ratio of ICON standard water, indexed as q^ICON. These tracers allow us to infer the relative importance of ocean and land evaporation – essentially the source of water in the model atmosphere – at all times. In addition to case studies, this is interesting because of the simplified implementation of isotopic processes during land evaporation in the current version of ICON-ART-Iso; see Sect. 2.3. The tracers of q^init provide information on the importance of the initialization at a certain time in the simulation.

3.1 An application of diagnostic water tracers: precipitation source regions

This section examines the moisture source regions of precipitation over ocean and land. Gimeno et al. (2012) give a review of the subject, while, e.g., Numaguti (1999), Van der Ent et al. (2010) and Risi et al. (2013) study this question by use of other models and in more detail.

Here, we use a decadal model integration. The simulation was initialized with ECMWF (European Centre for Medium-Range Weather Forecasts) Integrated Forecast System (IFS) operational analysis data on 1 January 2007, 00:00 UTC, to simulate 11 years on an R2B04 grid (≈160 km horizontal resolution). The time step was set to 240 s (convection called every second step) and output was saved on a regular $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid every 10 h in order to obtain values from all times of the day. Sea surface temperatures and sea ice cover were updated daily by linearly interpolating monthly data provided by the AMIP II project (Taylor et al., 2000). The first year is not considered as spin-up time of the model and the simulation is evaluated up to the end of 2017.

We look at the total precipitation P in Northern Hemisphere winter (December, January, February, denoted by DJF) and summer (June, July, August, denoted by JJA). Figure 1 displays zonal sums of P^init, P^ocn and P^lnd relative to standard water precipitation P^ICON as a function of latitude. The sum of precipitation that originates from convection is also given for each water species. The top panels give winter values, while the bottom panels display the results for summer months.

The area covered by ocean is not equally distributed over different latitudinal bands, which is the reason why ocean and land points are considered separately. Panels (b) and (e) show the fraction of precipitation that has fallen over the ocean relative to the total precipitation and the area fraction of the ocean in each latitudinal band. Despite the characteristics of the different seasons, which will be discussed in the following paragraphs, the latitudinal distribution of the ocean area fraction largely determines the overall fraction of rain that falls over the ocean or over land. This is why the other panels display values of P relative to the sum over each compartment, not to the total sum.

Figure 1 Fractional contributions of P^ICON, P^init, P^ocn and P^lnd to zonal sums of total precipitation for Northern Hemisphere winter (DJF, a–c) and summer (JJA, d–f) as a function of latitude. Outer panels show sums over the ocean and land grid points, respectively. Here, dashed lines indicate the contribution of convective precipitation for each source of atmospheric water. Center panels display the fraction of precipitation over the ocean relative to total precipitation (over land plus ocean) and the fraction of the area covered by ocean.

To evaluate the model simulation, we use values starting in the second year. At this time, the tropospheric moisture has been completely replaced by water that has evaporated during the model run. This is demonstrated by the values of P^init close to zero in all panels in Fig. 1. Technically, this means that the ternary solution of q^init, q^ocn and q^lnd that makes up q^ICON is practically reduced to a binary solution of only q^ocn and q^lnd. Other experiments show that this is already true after a few weeks (not shown).

During Northern Hemisphere winter over the ocean (Fig. 1a), the precipitation is strongly dominated by water that has evaporated from the ocean. Water from the land surface hardly reaches the ocean. Over land areas, the ocean is also the dominant source for precipitation, reaching more than 50 % at almost all latitudes. In the Northern Hemisphere midlatitude land areas, more than 70 % of the precipitated water originates from the ocean. The tropical and Southern Hemisphere land areas (in DJF) receive up to 40 % of precipitation from land evaporation. Most precipitation at tropical and subtropical latitudes over the ocean originates from convection (indicated by dashed lines). The role of convection is much smaller over land areas and again stronger in the Southern Hemisphere. Note that in a simulation with very high horizontal resolution (Klocke et al., 2017), more convective processes could have been directly resolved. In this specific case of a resolution close to 160 km, practically no convection is directly resolved by the model. It should therefore be considered that the amount of precipitation from convection only shows the importance of this parameterization in the simulations at this resolution.

The distribution of precipitation water sources is different in Northern Hemisphere summer (bottom row of Fig. 1). In summer, the Northern Hemisphere land areas (bottom right) supply themselves with a substantial fraction of the moisture that then precipitates. The importance of convection is increased in Northern Hemisphere summer with its maximum influence shifted into northern midlatitudes. Despite the larger moisture availability over the ocean, the far Northern Hemisphere land areas also supply the larger part of moisture that precipitates over the ocean in summer; see Fig. 1d.

These results are comparable to the studies by Numaguti (1999), Van der Ent et al. (2010) and Risi et al. (2013). While these studies look at regional differences, the latitudinal dependence is similar to the results presented here. This first application of ICON-ART-Iso – while no isotopologues are used – shows how diagnostic moisture tracers can be applied to better understand specific aspects of the atmospheric water cycle.

3.2 The multi-annual simulation compared to GNIP data

For a first validation of δD and δ¹⁸O values, we use the decadal ICON-ART-Iso model integration of the previous section and compare results to data taken from the GNIP network (IAEA/WMO, 2017; Terzer et al., 2013). In this section, we analyze δ values in total precipitation.

Five GNIP stations were chosen for their good data availability in the respective years, sampling different climate zones: Vienna in eastern Austria (48.2^∘ N, 16.3^∘ E) in central Europe, Ankara in central Anatolia (40.0^∘ N, 32.9^∘ E), Bangkok in tropical southern Asia (13.7^∘ N, 100.5^∘ E), Puerto Montt in central Chile (41.5^∘ S, 72.9^∘ W) and Halley station in Antarctica (75.6^∘ S, 20.6^∘ W). The closest grid point to each of these stations was taken from the model output and the multiyear mean of each calendar month was calculated for δD, δ¹⁸O and d-excess (d-excess $=\mathit{\delta }\text{D}-\mathrm{8}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) in precipitation, total precipitation P and 2 m temperature T_2m. The corresponding values are available from GNIP. Results are displayed in Fig. 2. The panels for total precipitation also include the fractional contribution to precipitation by ocean and land evaporation (see previous section). All panels (except for the precipitation amount) show the 1σ standard deviation range for model and measurement data.

Figure 2 Monthly mean data for five GNIP stations (left to right: Vienna, Ankara, Bangkok, Puerto Montt and Halley station). Variables listed from top to bottom: δD, δ¹⁸O, d-excess (δD−8 δ¹⁸O), total precipitation P_tot and 2 m temperature (T_2m). Plots showing P_tot also include the percentage of land and ocean evaporation in precipitation. The 1σ standard deviation interval is indicated by dashed lines (except P for readability).

For most stations, the seasonal cycle of precipitation is reproduced by the model. This includes the summer minimum for Ankara and the strong winter precipitation in Puerto Montt. Precipitation is underestimated for Bangkok, especially in Northern Hemisphere spring. For all stations, the influence of land evaporation is strongest in their respective summer. Vienna and Ankara show a decreasing influence of the ocean in winter, typical for a more continental climate. For Puerto Montt, located between the Pacific and the Andean mountain range, and for Bangkok, almost all precipitating water originates from the ocean.

The seasonal cycle of temperature is reproduced for all stations. Winter temperatures are too cold in this model configuration for all stations. This temperature bias can partly be explained by the fact that the altitude of all stations is higher in the model because of the coarse grid, e.g., 550 m for the grid point identified with Vienna versus 198 m for the GNIP station. Also, the measured temperatures are slightly higher than mean monthly ERA-Interim (Dee et al., 2011) 2 m temperatures for the corresponding grid points (not shown).

Despite some biases, the mean values of δD and δ¹⁸O are well reproduced by ICON-ART-Iso for all five stations. The seasonal cycle is captured correctly in the Northern as well as the Southern Hemisphere. Values of d-excess are also of similar magnitude. Model data are more variable than the measurements. However, the model data are mostly within the standard deviation range of measurements. This demonstrates the capability of ICON-ART-Iso to simulate climatological patterns. The seasonal cycle and regional differences in δD and δ¹⁸O are correctly reproduced in this climatological integration.

3.3 Comparison with IASI satellite data for a seasonal perspective

Here, we compare pairs of {H₂O,δD} retrieved from MetOp/IASI remote sensing measurements with data from two simulations. The section closely follows the case studies presented by Schneider et al. (2017), who compared IASI retrievals and data from the global hydrostatic model ECHAM5-wiso (Werner et al., 2011).

3.3.1 IASI satellite data and model post-processing

IASI (Infrared Atmospheric Sounding Interferometer) A and B are instruments onboard the MetOp-A and MetOp-B satellites (Schneider et al., 2016). They measure thermal infrared spectra in nadir view from which free-tropospheric {H₂O,δD} pair data are derived. As the satellites circle the earth in polar sun-synchronous orbit, each IASI instrument takes measurements twice a day at local morning (approximately 09:30) and evening (approximately 21:30). The measurements are most sensitive at a height of approximately 4.9 km. An IASI {H₂O,δD} pair retrieval method has been developed and validated in the framework of the project MUSICA (MUlti-platform remote Sensing of Isotopologues for investigating the Cycle of Atmospheric water). The MUSICA retrieval method is presented by Schneider and Hase (2011) and Wiegele et al. (2014) with updates given in Schneider et al. (2016).

Schneider et al. (2017) present guidelines for comparing model data to the remote sensing data. First, retrieval simulator software is used for simulating the MUSICA averaging kernel using the atmospheric state of the model atmosphere. The simulated kernel is then applied to the original model state (x) in order to calculate the state that would be reported by the satellite retrieval product ($\widehat{\mathit{x}}$, see Eq. 23).

$$\begin{array}{}\text{(23)}& \widehat{\mathit{x}}=\mathit{A}(\mathit{x}-{\mathit{x}}_{\text{a}})+{\mathit{x}}_{\text{a}}\end{array}$$

Here, A is the simulated averaging kernel and x_a the a priori state. The a priori value used in the retrieval process for 4.9 km is at $\mathit{\{}\mathrm{1780}\mathrm{ppm},-\mathrm{217.4}\mathrm{‰}\mathit{\}}$. This value represents the climatological state of the atmosphere. In the retrieval process, the satellite radiance measurements are used for estimating the deviation of the actual atmospheric state from the a priori assumed state, and it is important to note that the remote sensing retrieval product is not independent from the a priori assumptions (see Schneider et al., 2016, for more details). In Schneider et al. (2017), these guidelines have been followed for comparison of IASI data with ECHAM5-wiso model data. We use the same approach for comparisons to ICON-ART-Iso and our results can be directly compared to the results from the hydrostatic global model ECHAM5-wiso.

In order to compare ICON-ART-Iso measurements with IASI data, a simulation of 12 months is used, which was initialized on 5 November 2013. This simulation uses a finer resolution of R2B06, corresponding to roughly 40 km. Again, we use varying ocean surface temperatures and sea ice cover; see the specifications in Sect. 3.1. As in Schneider et al. (2017), two target time periods are investigated from 12–18 February and 12–18 August. As has been pointed out in Sect. 3.1, the amount of water remaining in the troposphere from initialization is negligible by using lead times of 3 months. For this study, model output was interpolated to a regular $\mathrm{0.36}{}^{\circ }\times \mathrm{0.36}{}^{\circ }$ grid, which is close to the 40 km resolution of the numerical ICON grid in the tropics. Output was written for every hour of simulation.

IASI observations are only available at cloud-free conditions. In order to exclude cloud-affected grid points in the ICON data, the total cloud cover simulated by ICON was used, denoted by C_clct. All points with C_clct>90 % were excluded. The parameter C_clct goes into saturation quickly and 90 % is reached even for thin clouds. Surface emissivity E_srf is a necessary input parameter for the retrieval simulator. In this first study, E_srf was set to 0.96 over land and 0.975 over the ocean. This is in accordance with the mean values as given by Seemann et al. (2008). In addition, Schneider et al. (2017) show in a sensitivity study that errors on the order of 10 % in this value have only a limited influence on the averaging kernels as simulated by the retrieval simulator. We follow the method outlined by Schneider et al. (2017) and use values only when the sensitivity metric s_err<0.05.

We examine results for different areas over ocean and over land using all data from the satellite and the model in the respective areas. The scatter of {H₂O,δD} is not shown directly. Instead, the figures show the isolines of relative normalized frequency, which is explained in Appendix A. In addition, Rayleigh fractionation curves are indicated in all figures. These are the same as those given by Schneider et al. (2017).

3.3.2 Seasonal and daily cycle

Seasonal and daily cycle are investigated in {H₂O,δD} space. The seasonal cycle is discussed for different regions over the central Pacific Ocean. The daily cycle is considered in the tropics and subtropics, also investigating differences between land and ocean areas.

Figure 3 Isolines of the relative normalized frequency distribution for pairs of δD and H₂O (see Appendix A for the method) after processing ICON-ART-Iso data with the IASI retrieval simulator of Schneider et al. (2017) (a, b) and IASI data for the same time (c, d). Data from morning overpasses are shown for 12 to 18 February (a, c) and 12 to 18 August (b, d) 2014 for different latitudinal bands over the Pacific Ocean (140^∘ E $<\mathit{\lambda }<\mathrm{220}{}^{\circ }$ E longitude). Contour lines are indicated at 0.2, 0.6 and 0.9 of the normalized distribution.

First, the seasonal cycle over the Pacific Ocean is examined by comparing the two target periods in different areas (140^∘ E $<\mathit{\lambda }<\mathrm{220}{}^{\circ }$ E longitude and different latitudinal bins). Results are presented in Fig. 3, which includes the exact latitudes. IASI data (bottom panels) show the specific characteristics of the different regions. H₂O content is highest for tropical air masses and lowest for the highest latitudes in February and August. At the same time, tropical air is least depleted in HDO, while the highest latitudes show the lowest values of δD, i.e., are more depleted. When comparing February and August values at each latitude, a clear seasonal signal appears everywhere except for the tropics: during summer of the corresponding hemisphere, the air is more humid and more depleted in HDO. The distributions seem to shift from season to season along a line perpendicular to those of the Rayleigh model. The distribution in the tropics shows a broadened shape in August.

The results of ICON-ART-Iso are shown in the top panels of Fig. 3. The latitudinal dependence is similar to IASI: high H₂O and δD in the tropics and lower values for midlatitudes. The range of values is also very similar. The seasonal cycle in H₂O and δD is also reproduced to some degree, especially in the subtropical latitudes. The most obvious differences to IASI results occur in the Northern Hemisphere midlatitudes in summer, which show less negative values of δD in the model than in the satellite data, especially for humid situations. In winter, this may also be the case, but there are only a few humid values simulated at all or available in the satellite dataset. In general, the model shows a similar behavior as ECHAM5-wiso, the results of which are presented by Schneider et al. (2017).

Figure 4 Isolines of the relative normalized frequency distribution for pairs of δD and H₂O (see Appendix A for the method) after processing ICON-ART-Iso data with the IASI retrieval simulator of Schneider et al. (2017) (a, b) and IASI data for the same time (c, d). (a, c) Data corresponding to morning and evening overpasses for the tropics (10^∘ S $<\mathit{\phi }<\mathrm{10}{}^{\circ }$ N, all longitudes, summer and winter simulation) over land and over the ocean. (b, d) Morning and evening overpasses for the subtropics (22.5^∘ S $<\mathit{\phi }<\mathrm{35}{}^{\circ }$ N, summer simulation) over land (Saharan desert region, 10^∘ W $<\mathit{\phi }<\mathrm{50}{}^{\circ }$ E) and Atlantic Ocean (50^∘ W $<\mathit{\phi }<\mathrm{30}{}^{\circ }$ W). Contour lines are indicated at 0.2, 0.6 and 0.9 of the normalized distribution.

For the daily cycle in the tropics and subtropics, land and ocean points are considered separately (Fig. 4; see caption for exact definition of the bins). IASI shows a clear signal of the daily cycle for both the tropics and subtropics over land (bottom panels of Fig. 4). There is no such signal over the ocean, where morning and evening distributions are almost identical. Over land, the water vapor in the tropics and subtropics is more depleted of HDO in the morning. There is also a daily cycle in H₂O in the tropics: during morning overpasses, H₂O values are higher than in the evening. Schneider et al. (2017) argue that this is due to the cloud filter, which removes areas of heavy convection in the evening. In the morning, the clouds have disappeared, but high humidity remains, especially in the lower troposphere. This may partly be due to evaporation of raindrops, which explains the enhanced depletion in HDO (Worden et al., 2007). Over the Sahara (the subtropical land area considered), the daily cycle is different: while mixing ratios of H₂O rise only slightly during the day, there is a strong increase in the HDO content in the evening. This behavior can be attributed to vertical mixing (Schneider et al., 2017).

The data retrieved from ICON-ART-Iso model simulations are shown in the top panels of Fig. 4. Tropical air (panel a) over the land shows slightly lower mixing ratios for H₂O than IASI. The humidity of tropical ocean points is better reproduced. The difference in δD is stronger for both areas, with δD values being too high in the model. There is no daily cycle in the tropics for ICON-ART-Iso. The subtropical mixing ratios (panel b) of H₂O over the ocean are similar to those in the tropics but cover a smaller range than those retrieved from IASI. The very humid and very dry parts of the IASI distribution are not reproduced by the model. δD values in ICON-ART-Iso are larger than in the IASI retrievals. As pointed out by Schneider et al. (2017), the daily cycle in IASI also manifests itself in the number of samples passing the IASI cloud filter and quality control. The IASI cloud filter removes much more evening observations than morning observations, meaning more cloud coverage in the evening than in the morning. In contrast, the ICON-ART-Iso cloud filter removes a similar amount of data for morning and evening; i.e., in the model morning and evening cloud coverage is rather similar. This may also influence the results.

To further analyze the influence of ocean and land areas, the analysis of the daily cycle is repeated, making use of the humidity tracers q^ocn and q^lnd. As has been pointed out in Sect. 3.1, q^init is negligible 3 months after initialization. To distinguish between grid points mostly influenced by ocean or land evaporation, we additionally use the following criteria to define ocean and land points: $q^{\text{ocn}}/q^{\text{ICON}}>\mathrm{0.9}$ for grid points over the ocean and $q^{\text{lnd}}/q^{\text{ICON}}>\mathrm{0.5}$ for land grid points predominantly affected by land evaporation. This investigation serves to showcase how the ocean and land evaporation tracers can be used and the threshold values are therefore arbitrary to some degree. The tracer fields of water evaporating from ocean and land have not been processed with the retrieval simulator, and instead values interpolated to 4.9 km are directly used.

Figure 5 As Fig. 4 for ICON-ART-Iso. In addition to the land–ocean mask, land data must pass the condition $q_{\mathrm{v}}^{\text{lnd}}/q_{\mathrm{v}}>\mathrm{0.5}$ and ocean data must pass $q_{\mathrm{v}}^{\text{ocn}}/q_{\mathrm{v}}>\mathrm{0.9}$.

The result is shown in Fig. 5 for the tropics and subtropics using the same method as for Fig. 4. The characteristics of the different regions show up much more clearly with the additional criteria. For the tropical ocean, the distribution of H₂O is similar, but the values are slightly more depleted in HDO. The distribution of pairs attributed to the land surface is reduced to values with relatively high humidity and enriched in HDO. The latter might be due to the signal of plant evapotranspiration, which is considered a non-fractionation process.

In the subtropics, the distributions over land change their shape completely and are separated from those over the ocean. The distribution for the subtropical ocean remains largely unchanged, becoming slightly more elongated with lower values in δD. For the land surface, the additional criterion strongly reduces the number of values that are considered. This implies that over the Saharan desert, air mostly influenced by land evaporation (50 % or more) is very dry and highly processed (low δD).

This section shows that ICON-ART-Iso is able to reproduce regional differences and the seasonal cycle of {H₂O,δD} of vapor in the lower troposphere. The additional water diagnostics are used to study the behavior of the model in more detail and will help investigate measured distributions in future studies.

3.4 Comparing with in situ IAGOS-CARIBIC measurements

In this section, we present a first case study, in which results of ICON-ART-Iso are compared to in situ measurements of δD taken by the IAGOS-CARIBIC passenger aircraft at 9–12 km of altitude. Two flights in September 2010 are considered, which took place a few days after the passage of the tropical cyclone Igor over the Atlantic Ocean. The full dataset of all δD measurements taken by IAGOS-CARIBIC in the tropics is also used as a reference.

3.4.1 IAGOS-CARIBIC data and model post-processing

In the European research infrastructure IAGOS-CARIBIC, a laboratory equipped with 15 instruments is deployed onboard a Lufthansa A340-600 for four intercontinental flights per month. Measurements of up to 100 trace gases and aerosol parameters are taken in situ and in air samples (Brenninkmeijer et al., 2007). δD is measured using the instrument ISOWAT (Dyroff et al., 2010). It is a tunable diode-laser absorption spectrometer that simultaneously measures HDO and H₂O at wave numbers near 3765 cm⁻¹ to derive δD in vapor. The instrument is calibrated based on regular measurements (each 30 min) of a water vapor standard with 500 ppm H₂O and $\mathit{\delta }\text{D}=-\mathrm{109}$ ‰. The δD offset is derived by considering the data of the driest 5 % of the air masses sampled during each flight, which is typically 4–8 ppm H₂O. At the flight altitude of 10–12 km, this is without exception lowermost stratospheric air (LMS), for which a δD of −600 ‰ is assumed (Pollock et al., 1980; Randel et al., 2012). An assumed uncertainty of this LMS value of 400 ‰ translates to a relevant uncertainty of 20 ‰ at 100 ppm H₂O. Due to further sources of measurement uncertainty, the data have a total flight-specific systematic uncertainty up to 100 ‰. The total uncertainty is humidity dependent, decreasing towards higher humidity (e.g., 100 ‰ at 80 ppm H₂O versus less than 20 ‰ at 500 ppm H₂O; see Christner, 2015, for more details).

The in situ IAGOS-CARIBIC data are suitable for the analysis of processes on small scales. δD measurements are available as 1 min means, which translates to a spatial scale of approximately 15 km. This horizontal resolution is finer than the chosen ICON-ART-Iso configuration (R2B06 corresponding to 40 km) and is therefore suitable for a case study validation. Unfortunately, the uncertainty of δD data at humidity below approximately 40 ppm H₂O is too high to be used for analysis. Because of the systematic total uncertainty (see above), we use mean δD values from two flights through similar conditions.

In this section, measurements from a return flight from Frankfurt to Caracas on 22 September 2010 are analyzed (IAGOS-CARIBIC flight nos. 309 and 310, taking off at 10:16 UTC in Frankfurt and 22:12 UTC in Caracas, respectively). The two flights crossed the Atlantic approximately 2 days after Hurricane Igor had passed the flight track. The storm caused large-scale lofting of tropospheric air masses and a moistening at flight level. The high humidity at flight level (9–12 km) allowed many accurate δD measurements to be taken.

An ICON-ART-Iso simulation was initialized with ECMWF IFS analysis data from 12 September 2010 and with the isotope values initialized as explained in Sect. 2.7. This corresponds to a 10-day forecast for the time of the two flights. In this case, not all tropospheric water from the initialization has been replaced by water evaporated during the simulation at the time of analysis. However, the δD values adjust to local meteorology within a few days, developing realistic horizontal and vertical gradients. The simulation was set up on an R2B06 grid (≈40 km) with a time step of 240 s (convection called every second step). The output with a frequency of one snapshot per 15 min was examined on a 0.5^∘ regular grid and interpolated linearly to the position of the aircraft. Figure 6 shows δD in water vapor in the upper troposphere roughly 24 h before the flights cross the Atlantic. The vortex signal of the hurricane clearly shows up in the δD field. The flight paths are also indicated in Fig. 6.

Figure 6 δD in water vapor on model level 54 (≈260 hPa) on 21 September, 21:00 UTC, the date prior to that of the IAGOS-CARIBIC flights. The storm is visible in the western half of the plotted area (center approximately at 30^∘ N, 50^∘ W). The flight paths of IAGOS-CARIBIC flights 309 and 310 are also indicated by two lines, for which the departure locations are emphasized.

3.4.2 Results for flights in tropical regions

In order to compare values influenced by Hurricane Igor, model and measurement data from flights 309 and 310 are considered in latitudes around the storm track only (0^∘ N $<\mathit{\phi }<\mathrm{30}{}^{\circ }$ N). For reference, all tropical δD values from the IAGOS-CARIBIC database are also examined. To create a comparable dataset from model data, the results of the decadal simulation presented in Sect. 3.1 were interpolated to the location of these measurements and are treated in the same manner. As in Sect. 3.3, the distribution of pairs of {H₂O,δD} are examined. The results are shown in Fig. 7.

The distribution of IAGOS-CARIBIC δD measurements is shown in the left panel of Fig. 7. The tropical measurement sample (blue contours) consists of all relevant measurements (23.5^∘ S $<\mathit{\phi }<\mathrm{23.5}{}^{\circ }$ N). While most tropical values are centered around −500 ‰ in δD and 100 to 150 ppm H₂O, there is also a tail towards more humid pairs in the distribution. The lower limit in δD follows the curves of Rayleigh fractionation. The measurement data from flights 309 and 310 (red contours) show different characteristics. The range in H₂O is similar to the maximum density of all tropical values, but the samples are more depleted in HDO. The humid branch is not continuous and the area of high humidity is sparsely populated. In general, both distributions are limited by the detection limit of 40 ppm in H₂O, while contour lines may reach slightly lower values because of the smoothing that is applied in processing the data (see Appendix A).

Model results are shown in the right panel of Fig. 7. For this figure, model data along the flight tracks are used only when accurate δD measurements are available. A limit of 40 ppm is also applied to the model data. The isolines stretching to lower value pairs again result from smoothing the data.

Figure 7 Isolines of the relative normalized frequency distribution (contours at 0.05, 0.4 and 0.9; see Appendix A for the method) of IAGOS-CARIBIC measurements (left) and ICON-ART-Iso model simulations for tropical samples, interpolated onto the paths of two IAGOS-CARIBIC flights 309–310 (right). The total number of data points in each distribution is given by the value N. Model data are considered only in locations with measurements, and the H₂O measurement limit of 40 ppm is also considered in model data.

The distribution from the tropical model sample is in some ways similar to the one by all tropical IAGOS-CARIBIC measurements (comparing the blue contours of the two panels): there is a tail towards high humidities and the upper limit of δD is roughly at −400 ‰, while the lower limit is given by the second Rayleigh curve. The model sample is 4 % lower in δD on average (H₂O reduced by 18 %). The mean for all values within the lowest density contour in {H₂O,δD} is at $\mathit{\{}\mathrm{193.5}\mathrm{ppm},-\mathrm{478.0}\mathrm{‰}\mathit{\}}$ for measurements, while it is at $\mathit{\{}\mathrm{175.4}\mathrm{ppm},-\mathrm{482.0}\mathrm{‰}\mathit{\}}$ for the model sample and $\mathit{\{}\mathrm{131.7}\mathrm{ppm},-\mathrm{481.6}\mathrm{‰}\mathit{\}}$ for measurements versus $\mathit{\{}\mathrm{72.7}\mathrm{ppm},-\mathrm{543.6}\mathrm{‰}\mathit{\}}$ in the case of the highest density contour.

The main characteristics of the distribution for the two flights following Hurricane Igor are captured by the model ICON-ART-Iso (comparing the red contours of the two panels) and the center of the distribution is below the full tropical sample. The model sample is slightly reduced in δD (1.1 %) when comparing to the full tropical sample, but 12.3 % more humid. From this simulation and these measurements alone, it is difficult to say if these discrepancies result from errors in the meteorological representation of the hurricane or in the physical parameterizations of the model. The good agreement between model and measurements in general is promising, while details will need to be examined in future studies.

By using the other three diagnostic moisture tracers (initialization water q^init and water evaporating from the ocean and land, q^ocn and q^lnd), the model results are examined further. The two transatlantic flights spent little time over land areas. Accordingly, $q_{\mathrm{v}}^{\text{lnd}}$ only reaches an average of 2.8 % for both flights, 19.9 % at maximum. An average of 47.2 % of the sampled water originates from the initialization (59.6 % maximum), while the remainder has evaporated from the ocean in the course of the simulation. Part of the discrepancies between model and measurements may thus result from the simplified representation of δD in the initial vapor field, while the influence of the approximated land surface isotope values remains limited.

This is analyzed in more detail in Fig. 8. Values of δD and H₂O along the flight paths are combined with information on the origin of the water that is sampled in the model. W stands for the ratio of vapor that originates from land or ocean evaporation or initialization, e.g., $W^{\text{init}}=q_{\mathrm{v}}^{\text{init}}/q_{\mathrm{v}}$. In Fig. 8, the scatter is color-coded by W^ocn. Because W^lnd is very low during most parts of the flights and especially so over the ocean, $W^{\text{ocn}}=\mathrm{1}-W^{\text{init}}$ is a good approximation. In Fig. 8, green indicates W^lnd>10 %, when the approximation of a binary solution is not valid.

Figure 8 Scatter of δD against H₂O from ICON-ART-Iso interpolated to the flight paths of IAGOS-CARIBIC flights 309 and 310. Color coding indicates the ratio of $W^{\text{init}}\approx \mathrm{1}-W^{\text{ocn}}$ in percent. Locations with W^lnd>10 % are marked in green. Values are considered for which p<280 hPa.

Figure 8 shows the strong influence of the ocean. More than 50 % of the sampled vapor originates from ocean evaporation for long parts of the flights. Values with the highest values of W^init (low W^ocn) mostly follow the course of the Rayleigh fractionation lines. The lowest and highest values in δD are reached when W^ocn is high (low values of W^init), indicating that the isotopologues in the model have seen many fractionating and transport processes. This includes air mass mixing, but also the microphysical and convective processes that are imprinted on the isotopologue ratio. The exact nature of these processes remains to be investigated in future studies.

4 Conclusions

We present ICON-ART-Iso, the isotope-enabled version of the global atmospheric model ICON (Zängl et al., 2015). We describe the model formulation as well as a set of evaluation studies. By using parts of the ICON-ART infrastructure (Schröter et al., 2018), the model is very flexible in terms of the simulated moisture tracers. These can be set to resemble either H₂O (tagged water) or the stable isotopologues HDO or ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ if fractionation is turned on. The physics of fractionation are largely based on the model COSMOiso (Pfahl et al., 2012). The first part of this article gives a detailed explanation of the parameterizations that have been implemented in ICON-ART-Iso to simulate the fractionation of water isotopologues.

We first evaluate tagged H₂O tracers of moisture evaporating from land and ocean to investigate the moisture sources of precipitation. This demonstrates the capabilities of ICON-ART-Iso to use tagged water as an additional diagnostic. The latitudinal dependence is similar to that given by other studies (Risi et al., 2013). The following three sections then investigated the performance of the model for the simulation of the isotopologues considering (i) multi-annual, (ii) regional and (iii) mesoscale applications.

For a multi-annual evaluation, the simulated isotopologues HDO and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ from a decadal simulation on a relatively coarse grid (160 km horizontal resolution) are compared to measurements taken from the network of GNIP stations (IAEA/WMO, 2017; Terzer et al., 2013). The model simulates δD and δ¹⁸O reasonably well, reproducing the seasonal cycle of δ¹⁸O and the range in d-excess for different stations in the Northern and Southern Hemisphere. This investigation presents a first climatological application.

Regional differences in pairs of {H₂O,δD} in lower free-tropospheric water vapor are then compared to data retrieved from IASI satellite measurements for a summer and winter case (Schneider and Hase, 2011; Schneider et al., 2016, 2017; Wiegele et al., 2014). The latitudinal dependence of these pairs is comparable to those from IASI retrievals. The seasonal cycle over the Pacific Ocean and the overall values are reproduced by the model in both seasons. The difference between land and ocean surfaces in the tropics and subtropics in the model is of similar magnitude as in the measurements. However, the daily cycle that is observed in the satellite data is not reproduced in the model. Overall, the performance is similar to that of ECHAM5-wiso (Schneider et al., 2017; Werner et al., 2011).

In a mesoscale application, a first comparison with in situ measurements uses δD in upper-tropospheric water vapor from two IAGOS-CARIBIC flights (Brenninkmeijer et al., 2007) transecting the Atlantic and from all tropical IAGOS-CARIBIC δD measurements (Dyroff et al., 2010). ICON-ART-Iso is able to reproduce the general features of the tropical IAGOS-CARIBIC dataset. The characteristics of the samples taken during two flights shortly after Hurricane Igor in September 2010 are also captured by the model.

In all three applications, the tagged evaporation water from ocean or land surfaces proves to be a valuable tool. It reveals a seasonal cycle in the precipitation water origin or shows the influence of initialization in the case of the comparison with IAGOS-CARIBIC data.

ICON-ART-Iso is a promising tool for future investigations of the atmospheric water cycle. This study demonstrates the flexibility of ICON-ART-Iso in terms of the setup for different diagnostics but also in terms of horizontal resolution and timescale. For future applications, it will be interesting to use a nudging of meteorological variables towards analysis data to facilitate comparisons with measurements in different case studies. Fractionation will be implemented in different microphysical schemes to make the model numerically more efficient and even better applicable for climatological questions. Due to its flexible setup, ICON-ART-Iso is ready to simulate tracers corresponding to ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{17}}\mathrm{O}$ or to be used as a test bed for new microphysical parameterizations.

Appendix A: Preparing the relative normalized frequency distributions

Section 3.3 and 3.4 show and discuss distributions of {H₂O,δD}. The scatter of {H₂O,δD} is not shown directly as the figures would be too cluttered. Instead, the normalized relative frequency is discussed, the isolines of which are shown in the different figures. This method has been adopted from Christner (2015). Figure A1 shows the scatter and the isolines of normalized relative frequency for the IAGOS-CARIBIC measurements of flights 309–310, which are discussed in Sect. 3.4.

To arrive at the isolines, the data are binned in H₂O and δD on a grid of 5 ppm ×5 ‰. In the case of IASI data, the data are binned in log₁₀H₂O (ppm) × δD on a grid of 0.05×5. Histogram counts are then interpolated onto a 1000×1000 grid. This is smoothed with a Gaussian filter with a standard deviation of 20 (15 in the case of IASI). These smoothed data are then normalized by the sum of all value pairs and normalized by the maximum value. Within this array of smoothed counts, isolines are drawn at 0.9, 0.4 and 0.05 (0.9, 0.6 and 0.2 in the case of IASI).

Figure A1 Scatter and isolines of the relative normalized frequency distribution for tropical (latitude $\mathit{\phi }<\mathrm{30}{}^{\circ }$ N) measurements of δD and H₂O from IAGOS-CARIBIC flights 309 and 310 (September 2010). The figure demonstrates how the isolines (indicated at 0.1, 0.4 and 0.9) relate to the underlying scatter.

Appendix B: Fractionation during evaporation of hydrometeors following Blossey et al. (2010)

For evaporation and equilibration of rain and melting hydrometeors, the parameterization following Stewart (1975) has been presented in this paper; see Sect. 2.5 and Eq. (6) therein. As an alternative, we have also implemented the parameterization by Blossey et al. (2010). For completeness, this parameterization is given in Eq. (B1). In order to make a comparison easy, Eq. (B2) again states the parameterization by Stewart (1975). As above, both are given in the formulation for the evaporation of rain. The notation is the same as in the main body text.

$$\begin{array}{ll}{\displaystyle }& {}^{\mathrm{h}}{S}_{\mathrm{r}}^{\text{evap, Blossey}}=A{\displaystyle \frac{{}^{\mathrm{h}}f}{{}^{\mathrm{l}}f}}{\displaystyle \frac{{}^{\mathrm{h}}D}{{}^{\mathrm{l}}D}}\\ \text{(B1)}& {\displaystyle }& {\displaystyle }\left[R_{\mathrm{r}}{\mathit{\alpha }}_{\text{liq}}^{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l},\mathrm{\infty }}^{*}-{}^{\mathrm{h}}{\mathit{\rho }}_{\mathrm{v}}\left(\mathrm{1}+B_{\mathrm{l}}\left(\mathrm{1}-{\mathit{\alpha }}_{\text{liq}}{\displaystyle \frac{R_{\mathrm{r}}}{R_{\mathrm{v}}}}\right)\right)\right]\text{(B2)}& {\displaystyle }& {}^{\mathrm{h}}{S}_{\mathrm{r}}^{\text{evap, Stewart}}=A{\left({\displaystyle \frac{{}^{\mathrm{h}}D}{{}^{\mathrm{l}}D}}\right)}^n\left[R_{\mathrm{r}}{\mathit{\alpha }}_{\text{liq}}^{\mathrm{l}}{\mathit{\rho }}_{\mathrm{l},\mathrm{\infty }}^{*}-{}^{\mathrm{h}}{\mathit{\rho }}_{\mathrm{v}}\right]\text{(B3)}& {\displaystyle }& {\displaystyle }A={\displaystyle \frac{\mathrm{4}\mathit{\pi }a{}^{\mathrm{l}}f{}^{\mathrm{l}}D}{\mathrm{1}+B_{\mathrm{l}}}}\text{(B4)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{l}}={\displaystyle \frac{{}^{\mathrm{l}}D{L}_{\mathrm{e}}^{\mathrm{2}}{}^{\mathrm{l}}{e}_{\mathrm{l},\mathrm{\infty }}^{*}}{k_{\mathrm{a}}{\mathcal{R}}_{\mathrm{v}}^{\mathrm{2}}{T}_{\mathrm{\infty }}^{\mathrm{3}}}}\end{array}$$

By the above formulation, the difference of the two parameterizations is easily accessible. While the empirical equation by Stewart (1975) introduces the exponent n, which is determined by measurements, the parameterization by Blossey et al. (2010) modulates the saturation difference by a factor determined by the actual isotopologue ratios in the hydrometeor and the surrounding vapor. The ratio of ventilation factors f is an additional tuning parameter in the parameterization by Blossey et al. (2010), which is set to 1 in the standard setup. A detailed comparison to the results of Stewart (1975) is postponed to a later study in order to keep this paper concise.