Nonequilibrium Fractionation During Ice Cloud Formation in iCAM5: Evaluating the Common Parameterization of Supersaturation as a Linear Function of Temperature

Supersaturation with respect to ice determines the strength of nonequilibrium fractionation during vapor deposition onto ice or snow and therefore influences the water isotopic composition of vapor and precipitation in cold environments. Historically, most general circulation models formed clouds through saturation adjustment and therefore prevented supersaturation. To match the observed isotopic content, especially the deuterium excess, of snow in polar regions, the saturation ratio with respect to ice (Si) was parameterized, usually by assuming a linear dependence of Si on temperature. The Community Atmosphere Model Version 5 (CAM5) no longer applies saturation adjustment for the ice phase and thus allows ice supersaturation. Here, we adapt the isotope‐enabled version of CAM5 to compute nonequilibrium fractionation in ice and mixed‐phase clouds based on Si from the CAM5 microphysics and use it to evaluate the common parameterization of Si. Our results show a wide range of Si predicted by the CAM5 microphysics and reflected in the simulated deuterium excess of Antarctic precipitation; this is overly simplified by the linear parameterization. Nevertheless, a linear function, when properly tuned, can reproduce the average observed relationship between δD and deuterium excess reasonably well. However, only the model‐predicted Si can capture changes in microphysical conditions under different climate states that are not due to changes in temperature. Furthermore, parametric sensitivity tests show that with the model‐predicted Si, water isotopes are more closely tied to the model microphysics and can therefore constrain uncertain microphysical parameters.


Introduction
Stable water isotopologues (H 16 2 O, HD 16 O, and H 18 2 O, hereafter referred to as stable water isotopes) have shown great potential as natural tracers of the global water cycle and have provided valuable insight into Earth's past and present climate (e.g., Dansgaard et al., 1993;Markle et al., 2017;Moyer et al., 1996). The basis for their wide range of applications is isotopic fractionation, which occurs during phase transitions and is caused by the different thermodynamic properties of the isotopes: heavier water isotopes have higher binding energies and therefore preferentially go to the compound in which the molecules are bound most [H]. Nonequilibrium fractionation occurs during phase transitions where diffusion is important, that is, (1) evaporation of water from the surface or from rain drops, when there is a strong gradient of humidity, and (2) ice and mixed-phase cloud formation, which usually occurs in an environment that is supersaturated with respect to ice (Rogers, 1979). Nonequilibrium fractionation is stronger for H 18 2 O than for HD 16 O and therefore leads to deviations from the ratio of 8:1, which are commonly quantified by the deuterium excess, d = D − 8 · 18 O (Dansgaard, 1964).
Due to the complexity of processes involving fractionation, measurements of stable water isotopes would be difficult to interpret without the help of numerical isotope models. Since Dansgaard (1964) developed a Rayleigh model simulating isotopic variations in an isolated air parcel to explain the temperature effect and the amount effect (i.e., increasing depletion of heavy isotopes with decreasing temperature and increasing precipitation amount), numerical models have been widely applied to link measured isotopic variations to physical processes. Following the pioneering work of Joussaume et al. (1984) stable water isotopes have been incorporated into many atmospheric general circulation models (AGCMs): iCAM (Lee et al., 2007;, ECHAMwiso (Hoffmann et al., 1998;Werner et al., 2011), GISS ModelE (Schmidt et al., 2005), HadCM3 , ICON-ART-Iso (Eckstein et al., 2018), IsoGSM (Yoshimura et al., 2008), and LMDZiso (Risi et al., 2010), among others. These AGCMs solve the full set of equations governing the dynamics and physics of the atmosphere on a global three-dimensional grid and simulate fractionation during all phase transitions in the hydrological cycle. In contrast to Rayleigh models, isotope-enabled AGCMs provide a spatially and temporally complete picture of global isotopic variations in water vapor and precipitation, limited only by their spatial resolution and the parameterization of unresolved processes. They have been widely applied to support interpretations of isotope measurements in paleoarchives. For example, Werner et al. (2000) showed with the help of ECHAMwiso that seasonality of precipitation explains the discrepancy between borehole and isotope-derived temperatures in Greenland, and Sime et al. (2009) showed with the help of HadCM3 that the relation between temperature and isotope ratios in Antarctic ice cores is nonlinear and that therefore interglacial climates may have been warmer than previously thought. At the same time, isotope tracers can be used to constrain physical parameterizations in the models that could otherwise not be detected because of compensating errors in the moisture, pressure, or temperature fields (Field et al., 2014;Risi et al., 2012). Historically, AGCMs use saturation adjustment to determine the rates of condensation, evaporation, deposition, and sublimation to or from both liquid and ice phase condensate. Thus, supersaturation is removed after each time step by converting excess vapor to liquid or ice clouds. However, no supersaturation means no nonequilibrium fractionation during ice cloud formation, and thus unrealistic deuterium excess values, especially at high latitudes. To solve this problem, supersaturation with respect to ice in these models is parameterized and accounted for in the microphysical transfers of the heavy water isotopes (HD 16 O, H 18 2 O) but not in the microphysical treatment of standard water (H 16 2 O). In most AGCMs, the same parameterization as the one used in Rayleigh models is adopted, which describes the saturation ratio with respect to ice (Si) as a linear function of temperature (T) in the form of Si = a + b · T (Jouzel & Merlivat, 1984;Petit et al., 1991). The intercept a is usually set to a value close to 1, implying saturation at T ≈ 0 • C, while the slope b is adjusted such that the simulated deuterium excess matches observations, usually of snow in Antarctica (e.g., Petit et al., 1991). The resulting values for b range from b = −0.002  to b = −0.007 (Schoenemann et al., 2014).
There are several limitations related to parameterizing supersaturation with respect to ice as a linear function of T: (1) A one-dimensional linear function cannot reflect a potential wide range of supersaturations at a given T; (2) since it is mostly tuned for present-day climate, the parameterization might not represent past climate conditions correctly, due to possible differences in the microphysics that influence supersaturation (e.g., different dust emissions or a different ratio of mixed-phase and ice clouds); and (3) since the isotope and microphysics parts of the model act independently in such a setup, an agreement between modeled

10.1029/2019MS001764
and observed isotope ratios reflects on the quality of the tuning process rather than the quality of the model (Mathieu et al., 2002).
The Community Atmosphere Model Version 5 (CAM5) (Neale et al., 2010) is one of a growing number of AGCMs whose microphysics scheme does not apply saturation adjustment for the ice phase and therefore allows supersaturation with respect to ice. This means that, although it is still used in the isotope-enabled version of CAM5 (iCAM5) , the parameterization of Si is not necessary. The aim of this study is therefore to adapt iCAM5 to compute nonequilibrium fractionation based on Si predicted by the microphysics scheme Morrison & Gettelman, 2008), and with the help of this new model version (1) evaluate the parameterization of Si as a linear function of T and (2) test the sensitivity of stable water isotopes to microphysical parameters influencing the model-predicted Si. Thereby, we focus on Antarctica, where supersaturation during ice and mixed-phase cloud formation is especially important for isotopes in precipitation and compare simulations in two different climates, present day and the last glacial maximum. Both iCAM5 versions are validated against observations of isotopes in surface snow and ice cores from Antarctica.

Model
CAM5 (Neale et al., 2010) is the atmospheric component of the National Center for Atmospheric Research Community Earth System Model (CESM) (Hurrell et al., 2013). For stratiform cloud physics it uses the macrophysics scheme by Park et al. (2014) and the two-moment microphysics scheme by Morrison and Gettelman (2008), modified to allow supersaturation with respect to ice by Gettelman et al. (2010). In pure ice clouds (T < −37 • C), ice crystals can form via heterogeneous immersion freezing on mineral dust or by the homogeneous freezing of sulfate (Liu & Penner, 2005;Liu et al., 2007). Ice crystal formation in mixed-phase clouds (−37 • C < T < 0 • C) occurs via heterogeneous deposition nucleation and condensation freezing (Meyers et al., 1992) and also through contact freezing (Young, 1974). Moist convection is separated into deep convection (Zhang & McFarlane, 1995) and shallow convection , and moist boundary layer turbulence is parameterized following Bretherton and Park (2009).
The implementation of stable water isotopes in CESM (Brady et al., 2019) follows the approach of previous modeling work (e.g., Joussaume et al., 1984;Hoffmann et al., 1998;Yoshimura et al., 2008). H 18 2 O and HD 16 O were added and tracked throughout the model's hydrological cycle. Both heavy isotopes experience equilibrium and nonequilibrium fractionation during phase transitions. The equilibrium fractionation factors ( eq ), which describe the ratios of saturation vapor pressures of the heavy isotopes (H 18 2 O, HD 16 O) and the light isotope (H 16 2 O), are parameterized as functions of T following Horita and Wesolowski (1994) for the liquid/vapor transition, and following Majoube (1971) and Merlivat and Nief (1967) for the ice/vapor transition. The ratios of diffusivities of the light and heavy isotopes (D∕D iso ), which are important for nonequilibrium fractionation, are taken from Merlivat (1978). In iCAM5 , fractionation during evaporation from the ocean is parameterized with the Craig and Gordon (1965) model using a wind speed-dependent formulation of the nonequilibrium fractionation factor (Merlivat & Jouzel, 1979). A similar approach is used for evapotranspiration from land in the isotope-enabled Community Land Model Version 4 (iCLM4)  that is coupled to iCAM5. Sublimation of surface snow or ice is assumed to occur without fractionation. The snow depth in iCLM4 is limited to a maximum snow water equivalent of H max = 1m, and any precipitation that would increase the snow depth to H > H max is routed directly to runoff (i.e., to the river component and ultimately to the ocean). This means that the isotopic composition of the snow pack is not updated in such regions. In the following, the implementation of fractionation during ice cloud formation, which is the focus of this study, is described in more detail. For a detailed description of other isotope parameterizations in iCAM5, see .

Fractionation During Ice Cloud Formation
Because of the low diffusivities of water molecules in ice and the short lifetime of ice crystals in the atmosphere, vapor deposition onto ice or snow is assumed to follow a Rayleigh process, which describes isotopic fractionation in an open system, where the condensate is immediately removed. The isotope ratio of the condensate (R c ) and the remaining vapor (R v ) are then given by where R v,0 is the initial isotope ratio in the vapor, f is the fraction of remaining vapor, and is the effective fractionation factor (including both equilibrium and nonequilibrium fractionation); is calculated following Jouzel and Merlivat (1984) and Blossey et al. (2010): where D and D iso are the diffusivities of the light and the heavy isotopes, respectively, Si is the saturation ratio with respect to ice, and eq is the equilibrium fractionation factor. In the default iCAM5 setup, Si is parameterized as a linear function of T, as in most previous modeling studies (e.g., Jouzel & Merlivat, 1984;Petit et al., 1991;Risi et al., 2010;Werner et al., 2011): with the tuning parameters a and b set to 1.0 and −0.002, respectively, in order to match observed precipitation deuterium excess over Antarctica .
To illustrate the impact of this parameterization, Figure 1 shows the evolution of D and deuterium excess in vapor and condensate described by the Rayleigh equations (equations (1) and (2) above) for three hypothetical air parcels with an initial vapor isotopic composition of D = −160‰ and 18 O = −20‰, and different parameterizations of Si as a linear function of T. The air parcels are assumed to cool from T = −10 • C to T = −50 • C with ∼10% of vapor condensing (to ice only) per 1 • C of cooling, such that the fraction of remaining vapor f = 2% at T = −50 • C.
We evaluate the effects of nonequilibrium fractionation in the D versus deuterium excess phase space, because it has been widely used in the past to analyze the effect of Si on isotopes (e.g., Petit et al., 1991;Risi et al., 2010;Werner et al., 2011). Since heavy isotopes condense more readily than light isotopes, the condensate has a higher D than the vapor it originates from in all three air parcels. As more vapor condenses, the removal of heavy isotopes leads to lower D in both vapor and newly forming condensate. The deuterium excess in vapor and condensate decreases first and then increases again. This curvature is due to an artifact of the traditional deuterium excess definition related to the nonlinearity of the scale (where the change of depends on itself) and could be avoided by using a definition based on the ln(R) scale (Dütsch et al., 2017;Markle et al., 2017). However, here we focus on the differences between the lines and therefore the curvature is of advantage, because it makes the differences more visible. Initially, at high D, a higher Si leads to a higher deuterium excess in the condensate and a lower deuterium excess in the remaining vapor, due to stronger nonequilibrium fractionation (i.e., HD 16 O diffuses more readily onto the ice crystals than H 18 2 O due to its larger diffusivity). However, because of the lower deuterium excess in the remaining vapor, the newly forming condensate has a lower deuterium excess as well, and after further deposition the Rayleigh lines of the condensate cross. For small fractions of remaining vapor, a higher Si therefore leads to a lower deuterium excess in the condensate. Such small fractions of remaining vapor are typically found at high latitudes or altitudes. Furthermore, a higher Si leads to an overall flatter slope between D and deuterium excess. Thus, Figure 1 demonstrates that Si has a large impact on the deuterium excess in both vapor and condensate, and therefore a correct representation of Si in isotope models is important.
In contrast to most AGCMs, the microphysics scheme of CAM5 allows Si > 100% , and the parameterization of Si is no longer necessary. Now vapor deposition onto ice and snow can be treated consistently for standard water (H 16 2 O) and the heavy isotopes (HD 16 O and H 18 2 O). In this study, we therefore adapt the computation of nonequilibrium fractionation in iCAM5 to use Si predicted by the microphysics scheme (hereafter referred to as Si real ) instead of the parameterized Si (hereafter referred to as Si param ).
Note that the default iCAM5 setup limits nonequilibrium fractionation during ice and mixed-phase cloud formation to T < −20 • C, whereas in reality it can occur at T < 0 • C if Si > 100%. For consistency this −20 • C limit is kept in the modified iCAM5 version but will be the focus of sensitivity tests (see section 2.2.2).

Table 1 Overview of Simulations
Control simulations (10 years) Sensitivity simulations (5 years) LGM, LGM, LGM, Note. Control simulations are run for present-day climate (PD) and the last glacial maximum (LGM). Sensitivity simulations are run for present-day climate only. T ini is the temperature threshold for nonequilibrium fractionation, naai is the number of active aerosols for ice nucleation, epsi is the Wegener-Bergeron-Findeisen time scale for the growth of ice crystals, and ai is the sedimentation velocity of ice crystals.

Simulations 2.2.1. Control Simulations
We run four control (ctrl) iCAM5 simulations, one for each of the Si versions (Si param and Si real ) used in two different climates (present day and last glacial maximum). All four simulations use the finite-volume dynamical core with 1.9 • N × 2.5 • E horizontal resolution and a hybrid sigma-pressure vertical coordinate system with 30 levels. iCAM5 is coupled to the isotope-enabled land surface model iCLM4  and the isotope-enabled sea ice model iCICE4 (Brady et al., 2019), with prescribed sea surface temperatures and sea ice concentrations. Present-day climate simulations use year 2000 CE climate forcings, and sea surface temperatures and sea ice concentrations from Hurrell et al. (2008). Last glacial maximum simulations use year 19050 BCE (21000 BP(1950)) climate forcings, and sea surface temperatures and sea ice concentrations from Zhu et al. (2017). Aerosol emissions are based off historical data from Lamarque et al. (2010) and Lamarque et al. (2011). Ocean isotope ratios are assumed to be constant in time and space ( 18 O = 0‰ and D = 0‰ in present-day climate, and 18 O = 1.13‰ and D = 9.10‰ during the last glacial maximum, which corresponds to the climatological global average from Zhu et al., 2017). For both climates, the simulations are branched from a 1-year spin-up simulation using Si param and run for 10 years (Table 1). In addition to the Si param simulations where the tuning parameter b in equation (4) is set to the default value (b = −0.002), we run two simulations for present-day climate and last glacial maximum with b = −0.004 and b = −0.006, respectively. The results of these simulations will be shown in Figure 6 and in the supporting information.

Sensitivity Tests
Since stable water isotopes are sensitive to Si, and Si is predicted by the microphysics scheme in the Si real simulations, isotopes can potentially be used to constrain microphysical parameters influencing Si in the model. This is especially useful because measurements of Si itself are sparse (Genthon et al., 2017). We select three microphysical parameters that significantly influence Si: the number of active aerosols for ice nucleation (naai), the Wegener-Bergeron-Findeisen time scale for the growth of ice crystals (epsi), and the sedimentation velocity of ice crystals (ai). To test the sensitivity of stable water isotopes to these parameters we run additional simulations for present-day climate, in which we scale the parameters by a factor , where = 0.5, 2.
Furthermore, we test the sensitivity of stable water isotopes to the temperature threshold (T ini ) below which nonequilibrium fractionation occurs for both Si versions, also in present-day climate. In addition to the control simulations using T ini = −20 • C, we run simulations with T ini = −10 • C and T ini = 0 • C. All sensitivity simulations are branched from the present-day 1-year spin-up simulation using Si param (and = 1, T ini = −20 • C) and run for 5 years (Table 1). black: the flux of standard water. For example, in the case of rain evaporation, rain is the source and vapor is the destination, or in the case of snow melting, snow is the source and rain is the destination.

T, Si, and Deuterium Excess in Vapor During Cloud Formation
Over Antarctica, most precipitation falls in the form of snow, which, in contrast to rain, does not experience strong fractionation during sublimation and thus carries its initial isotopic signal all the way to the ground. The deuterium excess in surface snow therefore strongly depends on the meteorological conditions (T and Si), as well as the isotopic composition of the vapor ( D v and 18 O v ) at cloud formation. To determine these variables, we add four new diagnostic water tracers to iCAM5, whose deposition and condensation fluxes to ice and liquid clouds in the microphysics scheme are equal to the fluxes of standard water multiplied by T, Si, D v and 18 O v , respectively ( Figure 2, red arrows). The ratio of tracer cloud mixing ratio to standard water cloud mixing ratio is thus equal to the value the variables had at the locations and times the cloud formed, weighted by the fluxes of standard water. By preserving this tracer ratio for fluxes between clouds and precipitation ( Figure 2, blue arrows), these values can be traced all the way to where precipitation falls. To ensure the tracer ratio corresponds to the conditions at cloud formation, the vapor phase of the tracers has to be equal to standard water vapor. This is done by setting tracer fluxes to and from vapor equal to the fluxes of standard water (Figure 2, black arrows). Note that with this method fluxes away from the source can differ from the fluxes to the destination (cf. Figures 2a and 2b), which violates mass conservation. For example, in the case of rain evaporation, the flux away from the source (rain) is equal to the flux of standard water multiplied by the tracer ratio in rain, while the flux to the destination (vapor) is equal to the flux of standard water. This is not a problem, because these tracers are purely diagnostic.

Observations
The simulated isotope ratios are compared with ice core measurements from five Antarctic sites:  Figure 3 shows the climatology of deuterium excess in precipitation for the four control simulations, and their differences (see supporting information Figures S1 and S2 for the Si param simulations with b = −0.004   In iCAM5, the deuterium excess during the last glacial maximum is generally lower over the ocean and higher over central Antarctica than in present-day climate (Figures 3c and 3f). The deuterium excess predicted by the different Si versions differs substantially over Antarctica (Figures 3g and 3h) and over Greenland (not shown), indicating the presence of nonequilibrium fractionation during ice and mixed-phase cloud formation and therefore the sensitivity of precipitation isotopes to Si in very cold regions. The model-predicted saturation ratio with respect to ice, Si real , leads to lower deuterium excess than Si param with b = −0.002, similar deuterium excess as Si param with b = −0.004, and higher deuterium excess than Si param with b = −0.006. The differences are most pronounced during the last glacial maximum. At the coast of Antarctica, the Si versions even show opposite trends between last glacial maximum and present-day climate. Si real produces lower deuterium excess values during the last glacial maximum than in present-day climate, while Si param produces higher or lower values, depending on the choice of b. The ice cores close to the coast of Antarctica also have lower deuterium excess values during the last glacial maximum than in present-day climate, in agreement with the Si real simulations. D and 18 O alone do not differ much between the Si param and Si real simulations relative to their spatial and temporal variations (see supporting information Figures S3 and S4). This is because they are primarily governed by equilibrium fractionation.  Figure 4 shows at which T and Si most Antarctic precipitation forms in iCAM5 (using the diagnostic water tracers described in section 2.3). Ninety percent of Antarctic precipitation forms at T between −42 and −14 • C in present-day climate ( Figure 4a) and between −50 and −20 • C during the last glacial maximum (Figure 4b). Si is mostly between 107% and 127% in present-day climate and between 107% and 130% during the last glacial maximum. Thus, despite much lower T, precipitation during the last glacial maximum does not form at significantly higher Si. The differences in occurrence frequency between the last glacial maximum and present-day climate also show a much more complex pattern than a simple shift to lower T and higher Si (Figure 4c).

Relationship Between T and Si
In both climates, Si is close to 100% at T ≈ 0 • C and increases with decreasing T as the ratio between the saturation vapor pressure over liquid water and ice increases. However, the increase of Si with decreasing T is not linear. The mean Si increases nonlinearly with decreasing T and is always higher than the default Si param with b = −0.002. Furthermore, for every depicted T there is a wide range of possible Si values, sometimes spanning more than 50%.
As Figure 5 shows, the differences between Si real and Si param also affect the deuterium excess. If Si real < Si param , the resulting deuterium excess in precipitation in the Si real simulations is always lower than in the Si param simulations with b = −0.002 (Figures 5a and 5d). If Si real > Si param , the signal is mixed: at relatively high T (> −37 • C), the resulting deuterium excess in precipitation is higher in the Si real simulations than in the Si param simulations with b = −0.002; at relatively low T (< −37 • C), it is generally lower in the Si real simulations than in the Si param simulations with b = −0.002. This pattern is the result of two additive effects: the difference in the strength of nonequilibrium fractionation, and the difference in the isotopic composition of the vapor from which clouds form. Nonequilibrium fractionation is stronger for higher Si and therefore the difference between deuterium excess in precipitation (d p ) and deuterium excess in vapor (d v ) is larger in the Si real simulations than in the Si param simulations where Si real > Si param , and smaller where Si real < Si param (Figures 5b and 5e). Because Si real is usually higher than Si param with b = −0.002, the stronger nonequilibrium fractionation removes more deuterium excess from the remaining vapor, and therefore, the deuterium excess in the remaining vapor is always lower in the Si real simulations than in the Si param simulations with b = −0.002 (Figures 5c and 5f). The opposite happens for Si param with b = −0.004 and b = −0.006 (see supporting information Figures S5 and S6), because in these cases Si real is usually lower than Si param . The effect of fractionation dominates at relatively high T (> −37 • C), leading to higher deuterium excess in the Si real simulations where Si real > Si param and vice versa, while the effect of the vapor dominates at relatively low T (< −37 • C), leading to lower deuterium excess in the Si real simulations compared to Si param with b = −0.002 for all Si (and higher deuterium excess compared to Si param with b = −0.004 and b = −0.006). Figure 6 shows the relationship between D and deuterium excess in precipitation in the control simulations and observations. In addition to the Si param simulation with the default b = −0.002, the simulations with b = −0.004 and b = −0.006 are shown in panels (c) and (d). Note the similarity between the polynomial fits and the Rayleigh lines in Figure 1 including the crossing points, indicating that the isotopic composition of precipitation over Antarctica can be approximated by a Rayleigh process.

Relationship Between D and Deuterium Excess
The simulated D and deuterium excess values interpolated to the ice core sites tend to be more enriched along the polynomial fits compared to the observations. This may be related to the relatively coarse resolution and consequently lower topography. The advantage of the D versus deuterium excess phase space is that we can still infer what the deuterium excess would be at the ice core sites if D was lower.
The simulations using Si real and Si param with b = −0.002 and b = −0.004, respectively, all produce a reasonable relation between D and deuterium excess in present-day climate (Figures 6a and 6c). In contrast, the simulation using Si param with b = −0.006 differs substantially from the observations, with too low deuterium excess values at lower D (< −300‰). All simulations underestimate deuterium excess at higher D (> −300‰), which might be related to evaporation from the ocean or land and is not further addressed in this study. The Si param simulation with b = −0.002 produces lower deuterium excess values at higher D and higher deuterium excess at lower D compared to the Si real simulation, which can be explained by the fact that the default parameterization mostly underestimates Si in iCAM5 as seen in Figure 4, leading to weaker nonequilibrium fractionation and a higher deuterium excess in the remaining water vapor (cf. Figure 1).
For the last glacial maximum there are only five observational data points, but they clearly fall within the range of values in the simulation using Si real (except for EDC), while they are overestimated by the simulation using Si param with b = −0.002 (Figures 6b and 6d). The simulation using Si param with b = −0.004 produces more reasonable D and deuterium excess values, but none of the Si param simulations capture the range of the ice core values as well as the Si real simulation does.
Thus, even though Figure 4 showed that Si is not a linear function of T, Si param can represent the 10-year average relationship between D and deuterium excess reasonably well. However, there is no guarantee that this is true in all situations. The large variability of Si at a given T might be important at shorter time scales, for example, for individual weather events, or in different regions, for example, in the upper troposphere. Therefore, we recommend using the more physically realistic Si real when possible.

Sensitivity Tests 3.4.1. Microphysical Parameters
This section addresses the sensitivity of the deuterium excess in Antarctic precipitation to the number of active aerosols for ice nucleation (naai), the Wegener-Bergeron-Findeisen time scale for the growth of ice crystals (epsi), and the sedimentation velocity of ice crystals (ai). By changing the properties of ice and mixed-phase clouds, these parameters directly influence Si. Furthermore, through feedbacks of clouds on climate, they indirectly influence T as well. In the Si param simulations, the deuterium excess is sensitive only to T, because Si is parameterized as a linear function of T. In the Si real simulations, the deuterium excess is sensitive to both T and Si. Figure 7 shows how T and Si are affected by the scaling of the microphysical parameters in present-day climate. Higher naai and epsi increase ice nucleation and deposition, respectively, and thus reduce supersaturation, leading to lower Si at all T (Figures 7a and 7b). Higher ai enhances the sedimentation of ice crystals, with a mixed effect on Si. Higher values of ai result in lower Si above −30 • C and higher Si below −30 • C (Figure 7c). On average, scaling all three microphysical parameters by = 2 compared to = 0.5 leads to higher T and lower Si. According to the Rayleigh lines in Figure 1, we therefore expect lower deuterium excess at high D and higher deuterium excess at low D in Antarctic precipitation for = 2 than for = 1 in both the Si param and the Si real simulations, and the opposite effect for = 0.5.
As Figure 8 shows, the lower Si and higher T in the simulations with = 2 indeed lead to a lower deuterium excess at high D and a higher deuterium excess at low D compared to the control simulation (vice versa for = 0.5). For all three parameters, the relation between D and deuterium excess is similar across the Si param simulations, because the effect of changes in on T alone is relatively small. In contrast, the deuterium excess in the Si real simulations is very sensitive to changes in Si caused by changes in . Since the initial isotopic composition of the vapor depends on processes that are not addressed in this study (e.g., evaporation from the ocean or land) we will not focus on the absolute values, but on the slope between D and deuterium excess instead. While the slope between deuterium excess and D is too steep compared to observations in the Si real control simulation, it is too flat in the Si real simulations with reduced naai and epsi, suggesting that the truth lies somewhere in between. Reducing ai leads to lower deuterium excess at low D, but with a similar slope between deuterium excess and D. Figure 9 shows how D and deuterium excess depend on the temperature threshold below which nonequilibrium fractionation occurs during ice and mixed-phase cloud formation (T ini ). Increasing T ini from −20 • C to −10 or 0 • C results in much lower deuterium excess values in the Si real simulations (Figure 9b), whereas the values in the Si param simulations do not depend strongly on T ini (Figure 9a). This can be explained by the difference between Si real and Si param at relatively high T (> −20 • C) (cf. Figure 4a). While Si param (with b = −0.002) only grows as large as 104% for T > −20 • C, Si real exceeds 110% in most instances when −20 • C < T < −10 • C. Therefore, Si real leads to stronger nonequilibrium fractionation than Si param , which leaves the remaining vapor (and consequently the newly forming condensate) depleted in deuterium excess. Interestingly, the effect of increasing T ini is highly nonlinear. The difference between the simulations with T ini = −20 • C and T ini = −10 • C is much larger than the difference between the simulations with T ini = −10 • C and T ini = 0 • C, suggesting that deposition onto ice and snow at T between −20 and −10 • C dominates the changes in both simulations with increased T ini .

Implications
Even though our results show that Si is clearly not a linear function of T, the Si param simulations, with the right tuning, produced a reasonable range of D and deuterium excess in Antarctic precipitation, both compared to observations and to Si real simulations. In part, this is not surprising, because the Si function had been tuned for the simulations to match the observations. Nevertheless, this means that models parameterizing Si as a linear function of T are still a valid tool for supporting the interpretation of isotope measurements in paleoarchives. However, there may be situations where the large Si variability independent of T is important, for example, for individual weather events. This variability can be represented only if Si is predicted by the microphysics scheme.
Furthermore, with Si predicted by the microphysics scheme, isotopes can serve as an additional observational constraint on microphysical parameters. Previous studies have shown that CAM5 tends to produce too much ice and too little supercooled liquid water in mixed-phase clouds (Cesana et al., 2015;Kay et al., 2016;Komurcu et al., 2014;Wall et al., 2017). This means that even if Si is simulated perfectly, the overestimated ice production can lead to too many nonequilibrium fractionation events and a lower deuterium excess in the remaining vapor, which is passed to precipitation forming further downstream. The ice bias in CAM5 can be caused by a too efficient Wegener-Bergeron-Findeisen process, that is, a too rapid growth of ice crystals at the expense of supercooled liquid water , and is often corrected for by decreasing the parameter epsi (e.g., Sagoo & Storelvmo, 2017;. Our simulation with reduced epsi has a flatter slope between D and deuterium excess, suggesting that epsi might indeed be too high in the default CAM5 setup. Other studies have shown that CAM5's ice nucleation scheme for mixed-phase clouds (Meyers et al., 1992) overestimates the concentration of ice nucleating particles at high latitudes (DeMott et al., 2010;Prenni et al., 2007;Xie et al., 2013), which also leads to too much ice in mixed-phase clouds. This is because the Meyers et al. (1992) scheme calculates ice nucleating particle concentration assuming a fixed dependence on T and Si based on measurements from the Sierra Nevada, where ice nucleating particles are much more abundant than at high latitudes, and does not take into account the spatial and temporal variability of ice nucleating particles. With most of Antarctic precipitation in our simulations forming in the mixed-phase cloud regime (−37 • C < T < 0 • C), naai is primarily predicted by the Meyers et al. (1992) scheme, and our results show that reducing naai improves the slope between D and deuterium excess as well. The overestimated ice fraction in mixed-phase clouds may also explain why allowing nonequilibrium fractionation at all T < 0 • C, instead of only at T < −20 • C, brings the simulated D and deuterium excess values much further away from observations: if the condensate was mainly liquid between T = −20 • C and T = 0 • C, this threshold would only have a minor impact.
With the implementation of a new ice nucleation scheme (Shi et al., 2015;Wang et al., 2014), as well as a new microphysics scheme (Gettelman, 2015), the ice bias has been greatly improved in newer versions of CAM (Bogenschutz et al., 2018). We therefore expect a better agreement between the modeled and observed D and deuterium excess without the need for reducing naai, epsi, or deactivating nonequilibrium fractionation at T > −20 • C in the isotope-enabled version of CAM6 that is currently under development.

Neglected Processes
One process that has been neglected in the discussion so far is evaporation from the ocean or land. Since evaporation is the only process involving strong nonequilibrium fractionation apart from ice and mixed-phase cloud formation, the deuterium excess is commonly used as a proxy for moisture source conditions (e.g., Jouzel et al., 1982;Uemura et al., 2012;Vimeux et al., 1999). From a Rayleigh perspective, moisture source conditions determine the initial isotopic composition of the air parcel, while meteorological conditions during cloud formation determine how the isotopic composition of the air parcel evolves. For Antarctic precipitation, moisture mainly originates from the ocean (Sodemann & Stohl, 2009). Previous studies have shown that the wind speed-dependent formulation by Merlivat and Jouzel (1979) that is used in iCAM5 for evaporation from the ocean does not represent nonequilibrium fractionation correctly (Bonne et al., 2019;Pfahl & Wernli, 2009;Uemura et al., 2010). This is why in this study we focus on the slope between D and deuterium excess and not on the absolute values. However, the isotopic composition of the initial vapor can also influence the slope, due to the nonlinearity of the scale (Dütsch et al., 2017;Markle et al., 2017). We test this with the Rayleigh equations (1) and (2) by adding −10‰ and +10‰ to the initial D of the air parcels, corresponding to an initial deuterium excess of −10‰ and +10‰, respectively. As Figure 10 shows, a higher initial deuterium excess in vapor unsurprisingly leads to a higher deuterium excess in the condensate, but the slope between D and deuterium excess is nearly independent of the initial values. We therefore expect very similar results with regards to the slope between D and deuterium excess for different formulations of fractionation during evaporation from the ocean or land.
By using the isotope ratios of precipitation from iCAM5, we also neglect potential postdepositional processes. Recent studies from Greenland and Antarctica show that isotopic exchanges between snow and water vapor (Steen-Larsen et al., 2014), fractionation during sublimation (Madsen et al., 2019;Ritter et al., 2016), and snow metamorphism (Casado et al., 2018) may alter the isotope ratios in surface snow, especially at low accumulation sites. Due to the snow depth limit of H max = 1m snow water equivalent in our version of CLM4 (see section 2.1), the simulated snow pack mainly reflects the initial isotopic composition instead of the signal from precipitation, and it is not possible to meaningfully account for postdepositional processes. In a revised version of CLM4 (van Kampenhout et al., 2017), the snow pack is allowed to refresh from the top, and any excess mass is removed from the lowest snow layer instead. Including this new treatment of snow in the isotope version of CLM4 (iCLM4)  will be the focus of future work.

Conclusions
The isotope version of the Community Atmosphere Model Version 5 (iCAM5) is one of a few isotope-enabled AGCMs whose microphysics scheme does not apply saturation adjustment for the ice phase and therefore allows supersaturation with respect to ice. In this study we adapted iCAM5 to compute nonequilibrium fractionation during ice and mixed-phase cloud formation based on the supersaturation predicted by the microphysics scheme (the real supersaturation, Si real ) instead of the commonly applied parameterization of supersaturation with respect to ice as a linear function of temperature (the parameterized supersaturation, Si param ).
A comparison between simulations using the real supersaturation and simulations using the parameterized supersaturation showed that a linear function oversimplifies the dependence of supersaturation with respect to ice on temperature and that differences between the real and the parameterized supersaturation are reflected in deuterium excess in Antarctic precipitation. The average relation between D and deuterium excess was nevertheless well reproduced by the simulations using the parameterized supersaturation, when properly tuned. Thus, a linear function of temperature tuned to match isotope measurements may be a reasonable approximation of supersaturation with respect to ice. However, for a more physically realistic treatment of nonequilibrium fractionation and to adequately represent microphysical changes that are independent of temperature, using the real supersaturation is preferable.
Our results also showed that with the real supersaturation stable water isotopes can constrain microphysical parameters that influence supersaturation with respect to ice in the model. Reducing the number of active aerosols for ice nucleation (naai) or the Wegener-Bergeron-Findeisen time scale for the growth of ice crystals (epsi) improved the relation between D and deuterium excess compared to observations, which is in agreement with previous findings showing that CAM5 overestimates the fraction of ice in mixed-phase clouds.
In summary, while the parameterization of supersaturation with respect to ice as a linear function of temperature may lead to reasonable results with the right tuning, using the real supersaturation to compute nonequilibrium fractionation ties the isotopes more closely to the model microphysics, which on the one hand facilitates interpretations of isotope measurements in paleoarchives, and on the other hand makes isotopes more useful observational constraints. Project P2EZP2_178439). We would like to acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Furthermore, we thank Hugh Morrison, Andrew Gettelman, Blaž Gasparini, and Bradley R. Markle for helpful discussions on microphysics and isotopes, and two anonymous reviewers for their constructive comments, which helped to improve the manuscript. iCESM is publicly accessible online (https:// github.com/NCAR/iCESM1.2). Model output data and scripts to create the figures are available online (https:// doi.org/10.5281/zenodo.3374014).