Controls on the water vapor isotopic composition near the surface of tropical oceans and role of boundary layer mixing processes

Understanding what controls the water vapor isotopic composition of the sub-cloud layer (SCL) over tropical oceans (δD₀) is a first step towards understanding the water vapor isotopic composition everywhere in the troposphere. We propose an analytical model to predict δD₀ motivated by the hypothesis that the altitude from which the free tropospheric air originates (z_orig) is an important factor: when the air mixing into the SCL is lower in altitude, it is generally moister, and thus it depletes the SCL more efficiently. We extend previous simple box models of the SCL by prescribing the shape of δD vertical profiles as a function of humidity profiles and by accounting for rain evaporation and horizontal advection effects. The model relies on the assumption that δD profiles are steeper than mixing lines, and that the SCL is at steady state, restricting its applications to timescales longer than daily. In the model, δD₀ is expressed as a function of z_orig, humidity and temperature profiles, surface conditions, a parameter describing the steepness of the δD vertical gradient, and a few parameters describing rain evaporation and horizontal advection effects. We show that δD₀ does not depend on the intensity of entrainment, in contrast to several previous studies that had hoped that δD₀ measurements could help estimate this quantity.

Based on an isotope-enabled general circulation model simulation, we show that δD₀ variations are mainly controlled by mid-tropospheric depletion and rain evaporation in ascending regions and by sea surface temperature and z_orig in subsiding regions. In turn, could δD₀ measurements help estimate z_orig and thus discriminate between different mixing processes? For such isotope-based estimates of z_orig to be useful, we would need a precision of a few hundred meters in deep convective regions and smaller than 20 m in stratocumulus regions. To reach this target, we would need daily measurements of δD in the mid-troposphere and accurate measurements of δD₀ (accuracy down to 0.1 ‰ in the case of stratocumulus clouds, which is currently difficult to obtain). We would also need information on the horizontal distribution of δD to account for horizontal advection effects, and full δD profiles to quantify the uncertainty associated with the assumed shape for δD profiles. Finally, rain evaporation is an issue in all regimes, even in stratocumulus clouds. Innovative techniques would need to be developed to quantify this effect from observations.

1 Introduction

1.1 What controls the water vapor isotopic composition?

The water vapor isotopic composition (e.g., $\mathit{\delta }D=\left(R/R_{\mathrm{SMOW}}-\mathrm{1}\right)\times \mathrm{1000}$ expressed in per mill, where R is the D∕H ratio and SMOW is the Standard Mean Ocean Water reference), has been shown to be sensitive to a wide range of atmospheric processes (Galewsky et al., 2016), such as continental recycling (Salati et al., 1979; Risi et al., 2013); unsaturated downdrafts (Risi et al., 2008, 2010a); rain evaporation (Worden et al., 2007; Field et al., 2010); the degree of organization of convection (Lawrence et al., 2004; Tremoy et al., 2014); the convective depth (Lacour et al., 2017b); the proportion of precipitation that occurs as convective or large-scale precipitation (Lee et al., 2009; Kurita, 2013; Aggarwal et al., 2016); vertical mixing in the lower troposphere (Benetti et al., 2015; Galewsky, 2018a, b), mid-troposphere (Risi et al., 2012b) or upper-troposphere (Galewsky and Samuels-Crow, 2014); convective detrainment (Moyer et al., 1996; Webster and Heymsfield, 2003); and ice microphysics (Bolot et al., 2013). It is therefore very challenging to quantitatively understand what controls the isotopic composition of water vapor.

A first step towards this goal is to understand what controls the water vapor isotopic composition in the sub-cloud layer (SCL) of tropical (30^∘ S–30^∘ N) oceans. Indeed, this water vapor is an important source moistening air masses traveling to land regions (Gimeno et al., 2010; Ent and Savenije, 2013) and towards higher latitudes (Ciais et al., 1995; Delaygue et al., 2000). It is also ultimately the only source of water vapor in the tropical free troposphere since water vapor in the free troposphere ultimately originates from convective detrainment (Sherwood, 1996), and convection ultimately feeds from the SCL air (Bony et al., 2008). Therefore, the water vapor isotopic composition in the SCL of tropical oceans serves as initial conditions to understand the isotopic composition in land waters and in the tropospheric water vapor everywhere on Earth. We focus here on the SCL because, by definition, there is no complication by cloud condensation processes.

The goal of this paper is thus to propose a simple analytical equation that allows us to understand and quantify the factors controlling the δD in the water vapor in the SCL of tropical oceans. So far, the most famous analytical equation for this purpose has been the closure equation developed by Merlivat and Jouzel (1979) (MJ79). This closure equation can be derived by assuming that all the water vapor in the SCL air originates from surface evaporation. The water balance of the SCL can be closed by assuming a mass export at the SCL top (e.g., by convective mass fluxes) and a totally dry entrainment into the SCL to compensate for this mass export. The MJ79 equation has proven very useful to capture the sensitivity of δD and second-order parameter d-excess to sea surface conditions (Merlivat and Jouzel, 1979; Ciais et al., 1995; Risi et al., 2010d). However, the δD calculated from this equation suffers from a high bias in tropical regions (Jouzel and Koster, 1996). This bias can be explained by the neglect of vertical mixing between the SCL and air entrained from the free troposphere (FT). The MJ79 equation can better reproduce surface water vapor observation when extended to take into account this mixing (Benetti et al., 2015, hereafter B15). This extension requires us to know the specific humidity (q) and water vapor δD of the entrained air. To get these values, they assume that the air entrained into the boundary layer comes from a constant altitude. However, this does not reflect the complexity of entrainment and mixing processes in marine boundary layers.

1.2 Entrainment and mixing mechanisms

Figure 1 summarizes our knowledge about these entrainment and mixing processes. In stratocumulus regions, clouds are thin and the inversion is just above the lifting condensation level (LCL). Air is entrained from the FT by cloud-top entrainment driven by radiative cooling or wind shear instabilities (Mellado, 2017), possibly amplified by evaporative cooling of droplets (Lozar and Mellado, 2015). Both direct numerical simulations (Mellado, 2017) and observations of tracers (Faloona et al., 2005) and cloud holes (Gerber et al., 2005) show that air is entrained from a thin layer above the inversion, thinner than 80 m and as small as 5 m. The boundary layer itself is animated by updrafts, downdrafts, and associated turbulent shells that bring air from the cloud layer downward (Brient et al., 2019; Davini et al., 2017).

In trade-wind cumulus regions, the cloudy layer is a bit deeper. Observational studies and large-eddy simulations have pointed out the important role of thin subsiding shells around cumulus clouds, driven by cloud-top radiative cooling, mixing, and evaporative cooling of droplets (Jonas, 1990; Rodts et al., 2003; Heus and Jonker, 2008; Heus et al., 2009; Park et al., 2016). This brings air from the cloudy layer to the SCL. Subsiding shells may also cover overshooting plumes of the cumulus clouds, entraining FT air into the cloud layer (Heus and Jonker, 2008).

In deep convective regions, unsaturated downdrafts driven by rain evaporation (Zipser, 1977) are known to contribute significantly to the energy budget of the SCL (Emanuel et al., 1994). Large-eddy simulations show that subsiding shells, similar to those documented in shallow convection, also exist around deep convective clouds (Glenn and Krueger, 2014). In the clear-sky environment between clouds, turbulent entrainment into the SCL may also play a significant role (Thayer-Calder and Randall, 2015).

Therefore, whatever the cloud regime, air entering the SCL from above may originate from either the cloud layer or the free troposphere, depending on the mixing mechanism. Therefore, in this paper in contrast with B15, we let the altitude from which the air originates, z_orig, be variable. We do not call it “entrained” air because entrainment sometimes refers to mixing processes through an interface (e.g.,  De Rooy et al., 2013; Davini et al., 2017), whereas air in the SCL may also enter through deep, coherent, and penetrative structures such as unsaturated downdrafts. We do not call it FT air either since it may originate from the cloudy layer.

Figure 1 Schematics showing the different types of clouds and mixing processes as a function of the large-scale circulation.

1.3 Goal of the article

To acknowledge the diversity and complexity of mixing mechanisms, we extend the B15 framework in several ways. First, we assume that we know the shape of δD profiles as a function of q. Second, we write the specific humidity of the air originating from above the SCL as a function of z_orig. Third, we account for rain evaporation and horizontal effects.

While B15 focused on observations during a field campaign, we also apply the extended equation to global outputs of an isotope-enabled general circulation model, with the aim to quantify the different factors controlling the δD variability in the tropics. The variable z_orig will emerge as an important factor. Therefore, we discuss the possibility that δD measurements at the near surface and through the lower FT could help estimate z_orig and thus the mixing processes between the SCL and the air above.

Note that we focus on δD only. Results for δ¹⁸O are similar. We do not aim at capturing the second-order parameter d-excess because our model requires some knowledge about free-tropospheric vertical profiles of isotopic composition. While δD is known to decrease with altitude (Ehhalt, 1974; Ehhalt et al., 2005; Sodemann et al., 2017), vertical profiles of d-excess are more diverse and less well understood (Sodemann et al., 2017). In addition, there is more need for an extension of MJ79 for δD than for d-excess since the effect of convective mixing is larger on δD than on d-excess (Risi et al., 2010d; Benetti et al., 2014).

2 Theoretical framework

2.1 Box model and budget equations

Building on Benetti et al. (2014) and B15, we consider a simple box representing the SCL (Fig. 2). We assume that the air comes from above (M) and from the incoming large-scale horizontal advection (F_adv) and is exported through the SCL top (N, e.g., turbulent mixing or convective mass flux) and by outgoing large-scale horizontal advection (F_adv,out). We assume that the SCL is at steady state. For example, its depth is constant. Since the SCL properties may exhibit a diurnal cycle (Duynkerke et al., 2004), this hypothesis restricts the application of this model to timescales longer than daily. The air mass budget of the SCL is thus

$$\begin{array}{}\text{(1)}& M+F_{\mathrm{adv}}=N+F_{\mathrm{adv},\mathrm{out}}.\end{array}$$

These fluxes also transport water vapor and isotopes. In addition, surface evaporation E and rain evaporation F_evap import water vapor and isotopes (Fig. 2).

Hereafter, to simplify equations, we use the isotopic ratio R instead of δD.

The SCL is usually well mixed (Betts and Ridgway, 1989; Stevens, 2006; De Roode et al., 2016). We thus assume that the humidity and isotopic properties are constant vertically and horizontally in the SCL. They are noted (q₀, R₀). The humidity and isotopic properties of the mass flux export N are thus also (q₀, R₀). The properties of the flux M are noted (q_orig, R_orig). The properties of the incoming air by horizontal advection are noted (q_adv, R_adv). For simplicity here we neglect the effect of horizontal gradients in humidity (i.e., q_adv=q₀), assuming that the main effect of horizontal advection on δD₀ arises from horizontal gradients in δD. Appendix C explains how R_adv can be calculated. At steady state, the water budget of the SCL is written

$$\begin{array}{}\text{(2)}& M\cdot q_{\mathrm{orig}}+E+F_{\mathrm{evap}}+F_{\mathrm{adv}}\cdot q_{\mathrm{0}}=(N+F_{\mathrm{adv},\mathrm{out}})\cdot q_{\mathrm{0}}.\end{array}$$

This model is consistent with SCL water budgets that have already been derived in previous studies (Bretherton et al., 1995), except that we consider steady state. This equation can be solved for q₀:

$$\begin{array}{}\text{(3)}& q_{\mathrm{0}}=q_{\mathrm{orig}}+{\displaystyle \frac{E+F_{\mathrm{evap}}}{M}}.\end{array}$$

The SCL humidity q₀ is thus sensitive to M, justifying that it can be used to estimate the mixing intensity or the “entrainment velocity” $w_{\mathrm{e}}=M/{\mathit{\rho }}_{\mathrm{0}}$ (ρ being the air volumic mass) (Bretherton et al., 1995).

At steady state, the water isotope budget of the SCL is written

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& M\cdot q_{\mathrm{orig}}\cdot R_{\mathrm{orig}}+E\cdot R_{\mathrm{E}}+F_{\mathrm{evap}}\cdot R_{\mathrm{evap}}+F_{\mathrm{adv}}\cdot q_{\mathrm{0}}\cdot R_{\mathrm{adv}}\\ & =(N+F_{\mathrm{adv},\mathrm{out}})\cdot q_{\mathrm{0}}\cdot R_{\mathrm{0}},\end{array}\end{array}$$

where R_E is the isotopic composition of the surface evaporation. It is assumed to follow the Craig and Gordon (1965) equation:

$$\begin{array}{}\text{(5)}& R_{\mathrm{E}}={\displaystyle \frac{R_{\mathrm{oce}}/{\mathit{\alpha }}_{\mathrm{eq}}-h_{\mathrm{0}}\cdot R_{\mathrm{0}}}{{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})}},\end{array}$$

where R_oce is the isotopic ratio in the surface ocean water, α_eq is the equilibrium fractionation calculated at the sea surface temperature (SST) (Majoube, 1971), α_K is the kinetic fractionation coefficient (MJ79), and h₀ is the relative humidity normalized at the SST ($h_{\mathrm{0}}=q_{\mathrm{0}}/q_{\mathrm{s}}(\mathrm{SST},P_{\mathrm{0}})$, where q_s is the saturation-specific humidity at SST and P₀ is the surface pressure).

We write the isotopic composition of the rain evaporation, R_evap, as

$$R_{\mathrm{evap}}={\mathit{\alpha }}_{\mathrm{evap}}\cdot R_{\mathrm{0}},$$

where α_evap is an effective fractionation coefficient. For example, if droplets are formed near the cloud base, some of them precipitate and evaporate totally into the SCL (e.g., in non-precipitating shallow cumulus clouds), then α_evap=α(T_cloud base). In contrast, if droplets are formed in deep convective updrafts after total condensation of the SCL vapor, and then a very small fraction of the rain is evaporated into a very dry SCL, then ${\mathit{\alpha }}_{\mathrm{evap}}=\mathrm{1}/\mathit{\alpha }\left(T_{\mathrm{SCL}}\right)/{\mathit{\alpha }}_{\mathrm{K}}$ (Stewart, 1975).

We note $\mathit{\eta }=F_{\mathrm{evap}}/E$ the ratio of water vapor coming from rain evaporation to that of surface evaporation, and $\mathit{\varphi }=F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}/E$ the ratio of water vapor coming from horizontal advection to that coming from surface evaporation. We note $\mathit{\beta }=R_{\mathrm{adv}}/R_{\mathrm{0}}$ the ratio of isotopic ratios of horizontal advection to that of the SCL.

Note that in all our equations, we assume that temperature and humidity profiles and all basic surface meteorological variables are known. We do not attempt to express either h₀ as a function of q₀ as in B15 or the q profile as a function of q₀. Our ultimate goal is to assess the added value of δD assuming that meteorological measurements are already routinely performed.

By combining all these equations, we get

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& R_{\mathrm{0}}=\\ & {\displaystyle \frac{(\mathrm{1}-r_{\mathrm{orig}})\cdot R_{\mathrm{oce}}/{\mathit{\alpha }}_{\mathrm{eq}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})\cdot r_{\mathrm{orig}}\cdot \left(\mathrm{1}+\mathit{\eta }\right)\cdot R_{\mathrm{orig}}}{(\mathrm{1}-r_{\mathrm{orig}})\cdot h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})\cdot \left(\mathrm{1}+\mathit{\eta }+(\mathrm{1}-r_{\mathrm{orig}})\cdot (\mathit{\varphi }\cdot (\mathrm{1}-\mathit{\beta })-\mathit{\eta }\cdot {\mathit{\alpha }}_{\mathrm{evap}}\right)}}\end{array},\end{array}$$

where $r_{\mathrm{orig}}=q_{\mathrm{orig}}/q_{\mathrm{0}}$ is the proportion of the water vapor in the SCL that originates from above.

An intriguing aspect of this equation is that the sensitivity to M disappears. In contrast to q₀, R₀ is not sensitive to M. Therefore, it appears illusory to promise that water vapor isotopic measurements could help constrain the entrainment velocity that many studies have striven to estimate (Nicholls and Turton, 1986; Khalsa, 1993; Wang and Albrecht, 1994; Bretherton et al., 1995; Faloona et al., 2005; Gerber et al., 2005, 2013). The lack of sensitivity of R₀ to M is explained physically by the fact that for a given q₀ and q_orig, if M increases, then E+F_evap increases in the same proportion to maintain the water balance. Therefore, the relative proportion of the water vapor originating from surface and rain evaporation to that coming from above, to which R₀ is sensitive, remains constant. Rather, since q and R vary with altitude, R₀ is sensitive to the altitude from which the air originates.

Figure 2 Schematics showing the simple box model on which the theoretical framework is based, and illustrating the main notations.

2.2 Closure if δD profile follows a Rayleigh distillation line

Equation (6) requires knowing q_orig and R_orig. B15 take these values from general circulation model (GCM) outputs at 700 hPa. In contrast, here we acknowledge the diversity and complexity of mixing mechanisms by keeping the possibility to take q_orig and R_orig at a variable altitude z_orig.

If the goal is to predict R₀ from z_orig, we can apply Eq. (6) if we know the q and δD vertical profiles. Conversely, if the goal is to predict z_orig from R₀, we can numerically solve Eq. (6) if we know the q and δD vertical profiles. No analytical solution exists in the general case, but a numerical solution can be searched for the z_orig based on Eq. (6). However, the existence and unicity of the solution is not warranted for all kinds of profiles (e.g., Appendix A).

In practice, full isotopic profiles are costly to measure. In addition, our goal is to develop an analytical model. Therefore, in the following we simplify the problem by assuming that the vertical profile of R follows a known relationship as a function of q. Measured vertical profiles of δD are usually bounded by two curves when plotted in a (q, δD) diagram (Sodemann et al., 2017): Rayleigh distillation curve and mixing line.

First, we explore the case of a Rayleigh distillation curve (Dansgaard, 1964), as in Galewsky and Rabanus (2016):

$$\begin{array}{}\text{(7)}& R_{\mathrm{orig}}=R_{\mathrm{0}}\cdot r_{\mathrm{orig}}^{{\mathit{\alpha }}_{\mathrm{eff}}-\mathrm{1}},\end{array}$$

where α_eff is an effective fractionation coefficient. Typically, q decreases with altitude, so R also decreases with altitude. However, in observations and models, vertical profiles of R can be very diverse (Bony et al., 2008; Sodemann et al., 2017). The water vapor may be more (Worden et al., 2007) or less (Sodemann et al., 2017) depleted than predicted by a Rayleigh curve using a realistic fractionation factor that depends on local temperature. Therefore, here we let α_eff be a free parameter larger than 1. Rather than assuming a true Rayleigh curve, we simply assume that R and q are logarithmically related. Effects of horizontal advection and rain evaporation on tropospheric profiles are encapsulated into α_eff.

Injecting Eq. (7) into Eq. (6), we get

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& R_{\mathrm{0}}={\displaystyle \frac{R_{\mathrm{oce}}}{{\mathit{\alpha }}_{\mathrm{eq}}}}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})\cdot \left((\mathrm{1}+\mathit{\eta })\cdot \frac{\mathrm{1}-r_{\mathrm{orig}}^{{\mathit{\alpha }}_{\mathrm{eff}}}}{\mathrm{1}-r_{\mathrm{orig}}}-\mathit{\eta }\cdot {\mathit{\alpha }}_{\mathrm{evap}}+\mathit{\varphi }\cdot (\mathrm{1}-\mathit{\beta })\right)}}\end{array}.\end{array}$$

A simpler form can be found if neglecting horizontal advection and rain evaporation effects ($\mathit{\varphi }=\mathit{\eta }=\mathrm{0}$):

$$\begin{array}{}\text{(9)}& R_{\mathrm{0}}={\displaystyle \frac{R_{\mathrm{oce}}}{{\mathit{\alpha }}_{\mathrm{eq}}}}\cdot {\displaystyle \frac{\mathrm{1}}{h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})\cdot \frac{\mathrm{1}-r_{\mathrm{orig}}^{{\mathit{\alpha }}_{\mathrm{eff}}}}{\mathrm{1}-r_{\mathrm{orig}}}}}.\end{array}$$

As a consistency check, in the limit case where the air coming from above is totally dry (r_orig=0), Eq. (9) becomes the MJ79 equation:

$$\begin{array}{}\text{(10)}& R_{\mathrm{0}}={\displaystyle \frac{R_{\mathrm{oce}}}{{\mathit{\alpha }}_{\mathrm{eq}}}}\cdot {\displaystyle \frac{\mathrm{1}}{h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})}}.\end{array}$$

Equation (8) tells us that whenever α_eff>1, R₀ decreases as r_orig increases (Fig. 3 red), i.e., as q_orig is moister. Therefore, R₀ decreases as z_orig is lower in altitude. This result may be counterintuitive, but can be physically interpreted as follows. If z_orig is high, mixing brings air with very depleted water vapor, but since the air is dry, the depleting effect is small. In contrast, if z_orig is low, mixing brings air with water vapor that is not very depleted, but since the air is moist, the depleting effect is large (Fig. 4a).

Figure 3 (red) shows that the range of possible δD values is restricted to −70 ‰ to −85 ‰. This explains why in quiescent conditions near the sea level in tropical ocean locations, the water vapor δD varies little (Benetti et al. (2014), Françoise Vimeux, personal communication, 2018). In the limit case where r_orig→1 (i.e., the air comes from the SCL top), $R_{\mathrm{0}}\to \frac{R_{\mathrm{oce}}}{{\mathit{\alpha }}_{\mathrm{eq}}}\cdot \frac{\mathrm{1}}{h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})\cdot {\mathit{\alpha }}_{\mathrm{eff}}}$ (L'Hôpital's rule was used to calculate this limit). This lower bound is not so depleted compared to the more depleted water vapor observed in regions of deep convection (e.g.,  Lawrence et al., 2002, 2004; Kurita, 2013). This is because when r_orig→1, the water vapor coming from above has a composition very close to that of the SCL, so the depleting effect is limited. In addition, surface evaporation strongly damps the depleting effect of mixing. Only rain evaporation or liquid–vapor exchanges (Lawrence et al., 2004; Worden et al., 2007) can further decrease R₀ (Appendix B).

Figure 3 (green) shows that the sensitivity to α_eff is relatively small but cannot be neglected. Therefore, predicting δD₀ requires having some knowledge about the steepness of the isotopic profiles in the FT. Rain evaporation and horizontal advection can have either an enriching or depleting effect, but do not qualitatively change the results (Fig. 3 purple and blue).

Now we consider the case of a mixing line. Detailed calculations in Appendix A show that the sensitivity to r_orig is lost. An infinity of FT end members can lead to the same δD₀ when mixed with the surface evaporation, as illustrated in Fig. 4b and analytically demonstrated in Appendix A. Our main results (more depleted δD₀ as r_orig increases, restricted range of δD₀ variations, relationship with z_orig) hold only for δD profiles that are steeper than a mixing line. This is the case for profiles that are intermediate between a Rayleigh and a mixing line, as is usually the case in nature (Sodemann et al., 2017) or in a general circulation model (Appendix D1).

Figure 3 δD₀ as a function of r_orig according to Eq. (9), with α_eff=α_eq as an example (red). For this illustrative purpose, we assume SST = 30 ^∘C, h₀=0.8, δD_oce=0 ‰, and $\mathit{\varphi }=\mathit{\eta }=\mathrm{0}$. The sensitivity to the effective fractionation factor α_eff (green) is shown. If rain evaporation is 25 % of surface evaporation (η=0.25), the solid pink and blue curves show the sensitivity to the effective fractionation factor α_evap. If the incoming water vapor by horizontal advection is 25 % of surface evaporation (ϕ=0.25), the dashed pink and blue curves show the sensitivity to the isotopic gradient quantified by β.

Figure 4 Idealized q−δD diagrams showing how the SCL water vapor δD is set, for the case where the tropospheric δD profile follows Rayleigh distillation (a) or a mixing line (b). For this illustrative purpose, we assume SST = 30 ^∘C, h₀=0.8, δD_oce=0 ‰, and $\mathit{\varphi }=\mathit{\eta }=\mathrm{0}$. In (a), the red curve shows the Rayleigh profile starting from the SCL and the purple curve shows the mixing line connecting the air coming from above to the surface evaporation, in the case $r_{\mathrm{orig}}=q_{\mathrm{orig}}/q_{\mathrm{0}}=\mathrm{0.7}$. The green curve shows the Rayleigh profile starting from the SCL and the blue curve shows the mixing line connecting the air coming from above to the surface evaporation, in the case $r_{\mathrm{orig}}=q_{\mathrm{orig}}/q_{\mathrm{0}}=\mathrm{0.5}$. One can see that when r_orig is lower, the mixing line is more curved, leading to more enriched values. In (b), the purple curve joins the SCL air and the air at all altitudes above the SCL. One can see that different values for r_orig can lead to the same value of δD in the SCL.

3 Model simulations, observations, and methods

3.1 LMDZ simulations

We use an isotope-enabled general circulation model (GCM) as a laboratory to test our hypotheses and investigate what controls the isotopic composition. We use the LMDZ5A version of LMDZ (Laboratoire de Météorologie Dynamique Zoom), which is the atmospheric component of the IPSL–CM5A coupled model (Dufresne et al., 2012) that took part in CMIP5 (Coupled Model Intercomparison Project; Taylor et al., 2012). This version is very close to LMDZ4 (Hourdin et al., 2006). Water isotopes are implemented the same way as in the predecessor LMDZ4 (Risi et al., 2010c). We use 4 years (2009–2012) of a simulation of the AMIP (Atmospheric Model Intercomparison Project) (Gates, 1992) that was initialized in 1977. The winds are nudged towards ERA-40 reanalyses (Uppala et al., 2005) to ensure a more realistic simulation. Such a simulation has already been described and extensively validated for isotopic variables in both precipitation and water vapor (Risi et al., 2010c, 2012a). The ocean surface water δD_oce is assumed constant and set to 4 ‰. The resolution is 2.5^∘ in latitude by 3.75^∘ in longitude, with 39 vertical levels. Over the ocean, the first layer extends up to 64 m, and a typical SCL extending up to 600 m is resolved by six layers. Around 2500 m, a typical altitude for the inversion for trade-wind cumulus clouds, the resolution is about 500 m.

For our calculations, we only use tropical grid boxes (30^∘ S–30^∘ N) over tropical oceans (>80 % ocean fraction in the grid box). In addition, to avoid numerical problems when estimating the effect of horizontal advection and rain evaporation, only grid boxes and days where E>0.5 mm d⁻¹ are considered. This represents 99.7 % of all tropical oceanic grid boxes.

Specific diagnostics for horizontal advection and rain evaporation are detailed in Appendix B and C.

3.2 STRASSE observations

We also apply our theoretical framework to observations during the STRASSE (Sub-Tropical Atlantic Surface Salinity Experiment) cruise that took place in the northern subtropical ocean in August and September 2012 (Benetti et al., 2014). This campaign accumulates several advantages that are important for our analysis: (1) continuous δD₀ measurements in the surface water vapor (17 m) at a high temporal frequency during 1 month (Benetti et al., 2014, 2015, 2017b), (2) associated surface meteorological measurements, including SST and h₀, (3) 22 radio soundings relatively well distributed over the campaign period and providing vertical profiles of altitude, temperature, relative humidity and pressure, (4) ocean surface water δD_oce measurements (Benetti et al., 2017a), (5) a variety of conditions ranging from quiescent weather to convective conditions, (6) on many vertical profiles, a well defined temperature inversion allows to calculate the inversion altitude.

We use δD₀ measurements on a 15 min time step. The measurements in ocean water were interpolated on the same time steps using a Gaussian filter with a width of 3 d. The radio-soundings are used together with all water vapor isotopic measurements that are within 30 min of the radio-sounding launch. Only profiles during the ascending phase of the balloon are considered because the descent phase is often located far away from the initial launch point (McGrath et al., 2006; Seidel et al., 2011).

3.3 Estimating the altitude from which the air originates

Here we explain how z_orig is estimated based on LMDZ outputs. First, we assume that the q and δD at 500 hPa (q_f, δD_f) belong to a Rayleigh distillation line starting from the surface with effective fractionation α_eff:

$${\mathit{\alpha }}_{\mathrm{eff}}=\mathrm{1}+{\displaystyle \frac{\ln (R_{\mathrm{f}}/R_{\mathrm{0}})}{\ln (q_{\mathrm{f}}/q_{\mathrm{0}})}}.$$

In a real field campaign, this assumption means that we do not need to measure the full vertical profile of δD, but only δD_f at a given free-tropospheric altitude (e.g., 500 hPa).

We checked that results are similar when defining the end member at 400 hPa rather than 500 hPa. However, the end member should be defined above 500 hPa to ensure that it is well above boundary layer processes. If the end member is defined below 500 hPa (e.g., 600 hPa), there are a few cases where q increases with altitude (q_f>q₀) due to horizontal advection or convective detrainment from nearby moister regions; meanwhile, δD decreases monotonically, leading to unrealistic values for α_eff.

Second, r_orig is estimated based on Eq. (9), using α_eff, α_eq, α_K , δD_oce, h₀, and δD₀ simulated by LMDZ.

Third, the altitude z_orig is estimated from r_orig. Using the q vertical profile, we find z_orig so that $q\left(z_{\mathrm{orig}}\right)=r_{\mathrm{orig}}\cdot q_{\mathrm{0}}$ (Fig. 5, red).

When estimating z_orig from observations, we follow the same methodology except that in absence of measurements for q_f and δD_f we assume a constant α_eff=1.07 based on LMDZ simulation and that α_eq, α_K , δD_oce, h₀, and δD₀ come from surface observations.

Note that r_orig and z_orig are not direct diagnostics from the simulation, but rather a posteriori estimates to match the simulated δD₀. Therefore, if assumptions underlying Eq. (9) are violated, then the estimate of r_orig, and subsequently z_orig, will be biased. The estimate of r_orig encapsulates the effect of mixing processes, but also all other processes that have been neglected in our theoretical framework, such as temporal variations in SCL depth, q₀, or δD₀ or vertical variations in q₀ or δD₀ within the SCL.

Figure 5 Schematics illustrating the typical structure of tropical marine boundary layers. The sub-cloud layer (SCL) extends from the surface to the lifting condensation level (LCL) and the cloud layer extends from the LCL to the inversion (z_i). EIS stands for the estimated inversion strength. Left: shape of the vertical profile in q (black) and q_s (green). Right: shape of the vertical profile in potential temperature θ, inspired by Wood and Bretherton (2006). The LCL, z_orig , $z_{\mathrm{orig},r_{\mathrm{orig}}=\mathrm{0.6}}$, and z_i altitudes defined in Sect. 3.4 are indicated.

3.4 Boundary layer structure diagnostics

Figure 5 illustrates the structure of a typical tropical marine boundary layer covered by stratocumulus or cumulus clouds (Betts and Ridgway, 1989; Wood, 2012; Wood and Bretherton, 2004; Neggers et al., 2006; Stevens, 2006). The cloud base corresponds to the lifting condensation level (LCL). Below is the well mixed SCL. Above is the cloud layer, topped by a temperature inversion. Above the inversion is the FT.

The LCL is calculated as the altitude at which the specific humidity near the surface equals the specific humidity at saturation of a parcel that is lifted following a dry adiabat (Fig. 5).

The temperature inversion is an abrupt increase in temperature that caps the boundary layer. Therefore, a method to automatically estimate its altitude z_i is to detect a maximum in the vertical gradient of potential temperature (Stull, 1988; Oke, 1988; Sorbjan, 1989; Garratt, 1994; Siebert et al., 2000). This method is sensitive to the resolution of vertical profiles (Siebert et al., 2000; Seidel et al., 2010). Therefore, we adapted this method in order to yield z_i values that best agree with what we would estimate from visual inspection of individual temperature profiles. In LMDZ, we calculate z_i as the first level at which the vertical potential temperature gradient exceeds 3 times the moist-adiabatic lapse rate. In observations, we calculate z_i as the first level at which the vertical potential temperature gradient exceeds 5 times the moist-adiabatic lapse rate because radio-soundings are noisier than simulated profiles.

Finally, we calculate z_orig(r_orig=0.6), which is the z_orig altitude if r_orig is set to 0.6. This usually coincides with the altitude of strong humidity decrease near the inversion (Fig. 5).

3.5 Averages and composites

All calculations are performed on daily values for LMDZ and on 15 min values for observations.

For LMDZ, when analyzing spatial and seasonal variability, seasonal averages are calculated at each grid box over tropical oceans by averaging all days of all years that belong to each season. Seasons are defined as boreal winter (December–January–February), spring (March–April–May), summer (June–July–August), and fall (September–October–November). For illustration purpose, all maps are plotted for boreal winter. Standard deviations are also calculated among all days of all years for each season.

The type of clouds and mixing processes depends strongly on the large-scale velocity at 500 hPa (ω₅₀₀, map shown in Fig. 6a), with shallow clouds in subsiding regions and deeper clouds in ascending regions (Fig. 1). Therefore, it is convenient to plot variables as composites as a function of ω₅₀₀ (Bony et al., 2004). To make such plots, we divide the ω₅₀₀ range from −30 to 50 hPa d⁻¹ into intervals of 5 hPa d⁻¹. In each given interval, we average all seasonal-mean values at all locations over tropical oceans for which seasonal-mean ω₅₀₀ belongs to this interval (e.g., Fig. 8a will be an example). Note that such composites are carried out on seasonal-mean ω₅₀₀ because cloud processes and their associated diabatic heating are tied to the large-scale circulation through energetic constraints (Yanai et al., 1973; Emanuel et al., 1994) that are best valid at longer timescales, otherwise, the energy storage term may become significant (e.g., Masunaga and Sumi, 2017). This is why ω₅₀₀ is generally averaged over a month or longer (e.g., Bony et al., 1997; Williams et al., 2003; Bony et al., 2004; Wyant et al., 2006; Bony et al., 2013). In addition, we primarily focus on understanding the seasonal and spatial distribution of δD₀.

The cloud cover strongly correlates with the inversion strength, which can be quantified by the estimated inversion strength (EIS;  Wood and Bretherton, 2006) (map shown in Fig. 6b) as a measure of inversion strength. We thus also plot variables as composites as a function of EIS. To make such plots, we divide the EIS range from −1 to 9 K into intervals of 0.5 K. In each given interval, we average all seasonal-mean values at all locations over tropical oceans for which seasonal-mean EIS belongs to this interval (e.g., Fig. 8b will be an example). Using seasonal-mean values is consistent with Wood and Bretherton (2006) and with the better link at longer timescales between cloud processes and the large-scale dynamical regime.

Figure 6 Maps of winter-mean ω₅₀₀ (a) and EIS (b) simulated by LMDZ.

3.6 Decomposition method for δD₀

To understand what controls the δD₀ spatiotemporal variations, δD₀ is decomposed into four contributions based on Eq. (8). First, we define $r_{\mathrm{orig},\mathrm{bas}}=\mathrm{0.3}$, ${\mathit{\alpha }}_{\mathrm{eff},\mathrm{bas}}=\mathrm{1.09}$, SST_bas=25 ^∘C, $h_{\mathrm{0},\mathrm{bas}}=\mathrm{0.7}$, ϕ_bas=0, η_bas=0, β_bas=1, and ${\mathit{\alpha }}_{\mathrm{evap},\mathrm{bas}}=\mathrm{1}$ as a basic state. We call $\mathit{\delta }{D}_{\mathrm{0},\mathrm{func}}(r_{\mathrm{orig}},{\mathit{\alpha }}_{\mathrm{eff}},\mathrm{SST},h_{\mathrm{0}},\mathit{\varphi },\mathit{\beta },\mathit{\eta },{\mathit{\alpha }}_{\mathrm{evap}})$ the function giving δD₀ as a function of r_orig, α_eff, SST, h₀, ϕ, β, η, and α_evap following Eq. (8), and $\mathit{\delta }{D}_{\mathrm{0},\mathrm{bas}}=\mathit{\delta }{D}_{\mathrm{0},\mathrm{func}}(r_{\mathrm{orig},\mathrm{bas}}$, α_eff,bas, SST_bas, h_0,bas, ϕ_bas, β_bas, η_bas, α_evap,bas). The relative contribution of r_orig to δD₀ is estimated as δD_0,func (r_orig, α_eff,bas, SST_bas, h_0,bas, ϕ_bas, β_bas, η_bas, ${\mathit{\alpha }}_{\mathrm{evap},\mathrm{bas}})-\mathit{\delta }{D}_{\mathrm{0},\mathrm{bas}}$. Similarly, the contributions of α_eff, SST, h₀, ϕ, and η to δD₀ are estimated as detailed in Table 1. All the contributions have the same units as δD₀ (‰). The sum of these components yields a quantity that is very close to the simulated δD₀, which confirms the validity of this linear decomposition. These components and their sum can be plotted as maps: Fig. 7 provides an example.

The relative contributions of each of these components to the δD variability are quantified by performing a linear regression of each of the components as a function of δD₀. If the correlation coefficient is significant for a given factor, then the slope quantifies the contribution of this factor to the variability of δD₀. The sum of all contributions may not always be 1 due to nonlinearity. Such a method has already been applied in previous studies (e.g., Risi et al., 2010b; Oueslati et al., 2016). The contributions to the seasonal spatial variability of δD₀ can be quantified by performing the regression among all locations and seasons. The contributions to the daily variability of δD₀ can be quantified by performing the regression among all days of a given season at a given location.

Table 1 Equations to calculate the relative contributions of r_orig, α_eff, SST, h₀, ϕ, and η to δD₀ and the physical meaning of these contributions.

3.7 Decomposition method for r_orig

To understand what controls r_orig, a similar method as for the decomposition of δD₀ can be applied. We can write r_orig as

$$\begin{array}{}\text{(11)}& r_{\mathrm{orig}}={\displaystyle \frac{h\left(z_{\mathrm{orig}}\right)\cdot q_{\mathrm{s}}\left(\stackrel{\mathrm{‾}}{T}\left(z_{\mathrm{orig}}\right)+\mathit{\delta }T\left(z_{\mathrm{orig}}\right),P\left(z_{\mathrm{orig}}\right)\right)}{q_{\mathrm{0}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{T}\left(z_{\mathrm{orig}}\right)+\mathit{\delta }T\left(z_{\mathrm{orig}}\right)=T\left(z_{\mathrm{orig}}\right)$ is the temperature at altitude z_orig, $\stackrel{\mathrm{‾}}{T}$ is the tropical-ocean-mean temperature profiles, h(z_orig) and P(z_orig) are the relative humidity and pressure at z_orig, and δT(z_orig) is the temperature perturbation compared to $\stackrel{\mathrm{‾}}{T}$. Therefore, the variability of r_orig is decomposed into the effect of four factors: q₀, z_orig, h(z_orig), and δT(z_orig). In practice, r_orig and z_orig are calculated following Sect. 3.3, and then Eq. (11) is applied.

4 Results from LMDZ

4.1 Decomposition of δD₀ variability

The spatial variations in δD₀ simulated by LMDZ (Fig. 7a) are characterized by depleted values near midlatitudes and in dry subsiding regions (e.g., off the coast of Peru and over other regions of oceanic upwelling) and regions of atmospheric deep convection (e.g., Maritime Continent). Consistently, δD₀ values exhibit a maximum for weakly ascending or subsiding regions: δD₀ decreases with increasing vertical velocity of both signs (Fig. 8a black); δD₀ decreases as EIS increases reflecting more stable, subsiding conditions (Fig. 8b black). This pattern is consistent with previous studies (e.g., Good et al., 2015). For the first time, we propose a theoretical framework to interpret this pattern, decomposing it into six contributions: r_orig, α_eff, SST, h₀, rain evaporation, and horizontal advection effects (Sect. 3.6). We check that the reconstructed δD₀ from the sum of its four contributions is very similar to the simulated δD₀ (Figs. 7b, 8 dashed black).

In ascending regions, the main contribution explaining the more depleted δD₀ in deep convective regions is that of α_eff (Figs. 7d, 8a red). α_eff is higher in more ascending regions (Fig. D1d). This means that the main factor depleting δD₀ in deep convective regions is the fact that the mid-troposphere is more depleted. This leads to a steeper gradient (higher α_eff), and thus a more efficient depletion by vertical mixing. This is consistent with deep convection depleting the water vapor most efficiently in the mid-troposphere (Bony et al., 2008). The second main contribution is that associated with r_orig (Figs. 7c, 8a green). r_orig is larger in deep convective regions (as explained in Sect. 4.2).

In subsidence regions, SST is the main factor controlling δD₀ (Figs. 7e, 8a pink): as subsidence is stronger, or as EIS increases, SST is colder, leading to larger α_eq and thus more depleted δD₀. Another important factor is h₀ (Figs. 7f, 8a purple): as subsidence is stronger, h₀ is drier, leading to more depleted δD₀. The contribution of r_orig is also a significant contribution to the depletion of δD₀ in the cold upwelling regions, for example off Peru or Namibia (Fig. 7c). The shallower boundary layer there is associated with higher r_orig.

The contribution of rain evaporation on δD₀ is minor compared to other contributions, except in the deepest convective regions (Fig. 7g). Rain evaporation is a slightly depleting effect in regions of strong deep convection and a slightly enriching effect in regions of moderate deep convection. When the fraction of raindrops that evaporate is small, isotopic fractionation favors evaporation of the lighter isotopologues. Therefore in convective, moist regions, rain evaporation has a depleting effect on the SCL (Worden et al., 2007). In contrast, in drier regions, rain evaporates almost totally. The evaporation flux thus has almost the same composition as the initial rain, which is more enriched than the water vapor.

The contribution of horizontal advection to δD₀ is significant only where isotopic gradients are the largest (Fig. C1h). Horizontal advection has slightly enriching in deep convective regions and depleting in coastal regions (e.g., off the coasts of California, Peru, Mauritania, Namibia, India, and Australia). For example, the Saharan layer off the northwestern African coast leads to a strong effect of horizontal advection (Lacour et al., 2017a).

From a quantitative point of view, we can decompose the δD₀ seasonal spatial variations into these different effects (Sect. 3.6). In regions of large-scale ascent, α_eff is the main factor explaining the δD₀ seasonal spatial variations (33 %), followed by rain evaporation (20 %) and r_orig (19 %; Table 2). In regions of large-scale descent, SST is the main factor explaining the seasonal spatial variations (54 %), followed by r_orig (29 %), h₀ (13 %), and α_eff (10 %) (Table 2). Note that the contribution of r_orig would be similar if we neglect rain evaporation and horizontal advection effects (Table 2).

The decomposition method can also be applied to decompose the δD₀ variability at the daily timescale at each location and for each season (Table 3). On average, in ascending regions, r_orig is the main factor (52 %), followed by rain evaporation (48 %) and α_eff (35 %). In subsiding regions, the effect of SST is muted due to its slow variability, and r_orig (82 %) becomes the main factor.

Overall, the results highlight the importance of r_orig as one of the main factors controlling the spatiotemporal variability of δD₀.

Figure 7 (a) Map of winter-mean δD₀ simulated by LMDZ. (b) Map of winter-mean δD₀ reconstructed as the sum of the four contributions. Tropical-mean δD₀ was added to compare with (a) on the same color scale. (c) Map of the contribution of r_orig on winter-mean δD₀ calculated from Eq. (9) (see Sect. 3.6). (d) Same as (b) but varies for α_eff. (e) Same as (b) but for SST. (f) Same as (b) but for h₀.

Figure 8 Composites as a function of ω₅₀₀ (a) and of EIS (b) of the seasonal averages of δD₀ simulated by LMDZ over all tropical ocean locations (black). Same for the sum of the contributions (black dashed and for each individual contribution to δD₀: r_orig varies (green), α_eff (dark red), SST (pink), h₀ (purple), rain evaporation (dashed brown), and horizontal advection (dashed blue). The tropical mean δD₀ was added to each contribution to plot on the same scale as simulated δD₀. The number of samples in each bin is indicated on a logarithmic scale on the right-hand side (dotted black line).

Table 2 Decomposition of the spatial and seasonal variation in δD₀ into its six contributions: effect of r_orig, α_eff, SST, h₀, rain evaporation, and horizontal advection. For each contribution, we show the correlation coefficient of the linear regression of the contribution as a function of δD₀. The analysis is performed separately for ascending and subsiding regimes. All seasons and locations over tropical oceans (30^∘ N–30^∘ S, ocean fraction >80 %, surface evaporation >0.5 mm d⁻¹) are considered. The threshold for the correlation coefficient to be statistically significant at 99 % is 0.15 or lower in all cases. We write correlation coefficient and slope values between brackets when they are not significant at 99 %.

Table 3 As in Table 4 but at the daily scale. The correlation coefficients and slopes are averaged over all seasons and locations over tropical oceans (30^∘ N–30^∘ S, ocean fraction >80 %), separately for ascending and subsiding regimes.

4.2 Decomposition of r_orig variability

Given the importance of r_orig in controlling the δD₀ variations, we now decompose r_orig into its four contributions: q₀, z_orig, h_orig, and δT_orig (Sect. 3.6). Spatially, r_orig is maximum in regions of strong large-scale ascent (Fig. 10a) such as the Maritime Continent (Fig. 9a) and in very stable regions (Fig. 10b) such as upwelling regions (Fig. 10a). We check that the reconstructed r_orig from the sum of its four contributions is very similar to the simulated r_orig (Figs. 9b, 10 dashed black).

In regions of strong large-scale ascent, r_orig is larger mainly because h_orig is larger (Figs. 9e, 10a pink). This is because the moister the FT, the higher the contribution of vapor coming from above to the vapor of the SCL, and thus the higher r_orig and the more depleted δD₀. This mechanism through which a moister FT leads to a more depleted δD₀ is consistent with that argued in B15. z_orig damps this effect: when convection is stronger and the FT moister, convection is also deeper, so the air originates from higher altitudes where the air is drier.

In very stable regions, r_orig is larger because q₀ is larger (Figs. 9c, 10b green), consistent with the drier conditions in these regions of large-scale descent. Note that this effect can be seen only in the most stable regions, but when considering all subsiding regions, the contribution is small (Table 2). r_orig is also larger because z_orig is lower in altitude (Figs. 9d, 10b red). As EIS increases, the boundary layers are shallower, the air comes from lower in altitude, r_orig is higher, and thus δD₀ is more depleted. This mechanism was not considered in B15 but our decomposition shows that it is a key mechanism driving r_orig and thus δD₀ variations in stable regions.

Quantitatively, in ascending regions, the main factor controlling the seasonal spatial variations in r_orig is h_orig (182 %), dampened by z_orig (−67 %) (Table 4). In descending regions, the main factor is also h_orig (96 %), followed by z_orig (41 %) (Table 4). At the daily scale, the same two factors dominate the variability of r_orig: h_orig and z_orig contribute to 78 % and 39 % of r_orig variations on average over ascending regions and to 118 % and 39 % on average over descending regions (Table 5).

Figure 9 (a) Map of winter-mean r_orig simulated by LMDZ. (b) Map of winter-mean r_orig reconstructed as the sum of the four contributions. Tropical-mean r_orig was added to compare with (a) with the same color scale. (c) Map of winter-mean r_orig calculated from Eq. (11) if only q₀ varies (see Sect. 3.6). (d) Same as (b) but if only z_orig varies. (e) Same as (b) but if only h(z_orig) varies. (f) Same as (b) but if only δT_orig varies.

Figure 10 Composites as a function of ω₅₀₀ (a) and of EIS (b) of the seasonal averages of r_orig simulated by LMDZ over all tropical ocean locations (black). Same for the sum of the contributions (black dashed) and for each individual contribution to r_orig:q₀ varies (green), z_orig (dark red), h(z_orig) (pink), and δT_orig (purple). The number of samples in each bin is indicated on a logarithmic scale on the right-hand side as bars.

Table 4 Decomposition of the spatial-seasonal variation in r_orig into its four contributions: effect of q₀, z_orig, h_orig, and δT_orig variations. For each contribution, we show the correlation coefficient of the linear regression of the contribution as a function of r_orig. The analysis is performed separately for ascending and subsiding regimes. All seasons and locations over tropical oceans (30^∘ N–30^∘ S, ocean fraction >80 %) are considered. The threshold for the correlation coefficient to be statistically significant is 0.15 or lower in all cases. We write correlation coefficient and slope values between brackets when they are not significant at 99 %.

Table 5 As in Table 4 but at the daily scale. The correlation coefficients and slopes are averaged over all seasons and locations over tropical oceans (30^∘ N–30^∘ S, ocean fraction >80 %), separately for ascending and subsiding regimes.

4.3 Estimating altitude z_orig

Estimated altitude z_orig is at a minimum in dry subsiding regions, especially in upwelling regions (Figs. 11a, and 12), corresponding to regions with the strongest inversion (Fig. 11). This contributes to the depleted δD₀ in these regions.

As explained in Sect. 3.3, our estimate of z_orig may be artificially biased due to the neglect of some processes in our theoretical framework. Ideally, to check whether z_orig really physically represents the altitude from which the air originates, additional model experiments where water vapor from different levels are tagged (Risi et al., 2010b) would be needed. While we leave this for future work, we check whether z_orig estimates are consistent with what we expect based on what we know about mixing processes in the marine boundary layers. We expect that in stratocumulus regions, air originates from a very shallow (a few tens of meters) layer above the inversion, whereas the mixing processes may be more diverse, and possibly deeper in the FT, as the boundary layer deepens (Fig. 1).

To check whether estimated z_orig is consistent with this picture, we compare z_orig to z_orig(r_orig=0.6) (z_orig that we would estimate is r_orig was set constant to 0.6) and z_i (Sect. 3.4), which are measures of the altitude of the humidity drop and temperature inversion, respectively. As expected from Fig. 1, they are minimum in dry upwelling regions, intermediate in trade-wind regions, and maximum values in convective regions (Figs. 11c–d, 12 green, blue). Therefore, the low z_orig in upwelling regions reflects the low z_i. Consistently, in subsiding regions, z_orig correlates well with $z_{\mathrm{orig},r_{\mathrm{orig}}=\mathrm{0.6}}$ (correlation coefficient of 0.52, statistically significant beyond 99 %). If we focus on very stable regions only (EIS >7 K), z_orig correlates well with both z_orig(r_orig=0.6) and z_i (correlation coefficient of 0.58 and 0.52, respectively, statistically significant beyond 99 %). The altitude z_orig is a few meters above the inversion in stratocumulus regions, and up to 1 km above the inversion in cumulus and deep convective regions (Fig. 12), consistent with our expectations from Fig. 1. This lends support to the fact that at least in subsiding regions, our isotope-based z_orig estimate effectively reflects the origin of air coming from above.

In ascending regions, in contrast, z_orig does not correlate significantly with z_orig(r_orig=0.6) or z_i. This may indicate either that our z_orig estimate is biased by neglected processes such as rain evaporation or that in deep convective regions the origin of FT air into the SCL is very diverse due to the variety of mixing processes (Fig. 1).

Figure 11 (a) Map of winter-mean z_orig estimated from δD₀ simulated by LMDZ. (b) Same as (a) but z_orig that we would estimate if r_orig were constant set to 0.6, z_orig(r_orig=0.6). (c) Same as (a) but for z_i simulated from LMDZ. Only days when EIS >2 K are considered, otherwise z_i is difficult to estimate. (d) Same as (a) but for LCL simulated by LMDZ.

Figure 12 Composites as a function of EIS of seasonal mean of z_orig (black), z_orig(r_orig=0.6) (green), z_i (blue), and LCL (red). The composite profiles of cloud cover are also shown, showing deep clouds when EIS is close to 0 and the shallowest clouds when EIS is largest.

5 Results from observations

To check whether our results obtained with LMDZ are realistic, we apply our methods to the measurements gathered during the STRASSE campaign. For simplicity and in absence of all necessary measurements, here we neglect the effects of rain evaporation and horizontal advection.

Throughout the cruise, δD₀ shows a large variability, ranging from around −75 ‰ in quiescent conditions to −120 ‰ during the two convective conditions (Benetti et al., 2014) (Fig. 13a red). Variability in r_orig is the major factor contributing to this variability (58 %) (Fig. 13a green, Table 6). This crucial importance of mixing processes is consistent with B15.

During the two convective events, the estimated r_orig saturates at 1 (Fig. 13b). This proves that r_orig estimated in these conditions is biased high because it encapsulates the effect of neglected processes, i.e., depletion by rain evaporation. Equation (9) is not valid in this case. In addition, at the scale of a few hours, the steady-state assumptions may be violated. Rain evaporation may strongly deplete the SCL before surface evaporation has the time to play its dampening role, hence the possibility of reaching very low δD₀ that cannot be predicted even when considering rain evaporation (Appendix B).

During the rest of the cruise, the main factors controlling the r_orig variability are z_orig (90 %) and h_orig (70 %). The importance of FT humidity in controlling r_orig was already highlighted in B15. However, in their paper, the variability in z_orig was neglected, whereas it appears here as the main factor.

Through September, the cruise goes from a shallow boundary layer in early September to deeper boundary layers with higher inversions, before reaching the convective conditions (Fig. 13c). Consistently with this deepening boundary layer, the air originates from increasingly higher altitudes. Remarkably, there are 6 d when z_orig coincides with z_i with a root-mean-square error of 31 ‰ and correlation coefficient of 0.996 (Fig. 13c). This indicates that the air comes exactly from the inversion layer. When recalling that z_orig and z_i are estimated from completely independent observations, the coincidence is remarkable and lends support to the fact that on these days, our z_orig estimate is physical. However, there remain 9 d when z_orig is much higher than z_i. This may reflect more penetrative downdrafts as we approach deeper convective regimes. But it may also be an artifact of our neglect of horizontal advection. For example, on these days which are characterized by lower h₀, neglecting the advection of enriched water vapor from nearby regions with higher h₀ could be misinterpreted as lower r_orig and thus higher z_orig.

Figure 13 (a) Time series of δD₀ observed during the STRASSE cruise, together with its four contributions. The δD of the surface ocean water is also plotted with the scale on the right. (b) Time series of r_orig estimated from observations during the STRASSE campaign, together with its four contributions. (c) Times series of z_orig, z_orig(r_orig=0.6), LCL, and z_i estimated from the STRASSE observations.

Table 6 Same as Table 2 but for the STRASSE observations. Linear regressions are calculated among 1977 data points.

Table 7 Same as Table 4 but for the STRASSE observations. Linear regressions are calculated among 55 data points, so that correlation coefficients above 0.35 are statistically significant at 99 %.

6 Discussion: what can we learn from water isotopes on mixing processes?

We have shown in the previous section that one of the main factors controlling δD₀ at the seasonal spatial and daily scales is the proportion of the water vapor in the SCL that originates from above (r_orig) and that one of the main factors controlling r_orig is the altitude from which the air originates (z_orig). In turn, could we use water vapor isotopic measurements to constrain z_orig? This would open the door to discriminating between different mixing processes at play (Fig. 1). Since mixing processes are crucial to determine the sensitivity of cloud fraction to SST (Sherwood et al., 2014; Bretherton, 2015; Vial et al., 2016), such a prospect would allow us to improve our knowledge of cloud feedbacks, and hence of climate sensitivity.

With this in mind, we assess the errors associated with z_orig estimates from δD₀ measurements, and discuss whether they are small enough for z_orig estimates to be useful. In stratocumulus clouds where the air is believed to originate from the first few tens of meters above cloud top (Faloona et al., 2005; Mellado, 2017), z_orig estimates are not useful if the errors are larger than a few tens of meters, e.g., 20 m. In cumulus clouds where mixing processes are more diverse and possibly deeper (Fig. 1), z_orig estimates may be useful if errors are of the order of 80 m.

Let us assume that we have a field campaign where we measure δD₀, surface meteorological variables, temperature and humidity profiles (e.g., radio soundings), and a few δD profiles (e.g., by aircraft). This is what we can expect for example from the future EUREC4A (Elucidating the role of clouds-circulation coupling in climate) campaign to study trade-wind cumulus clouds (Bony et al., 2017). Below we quantify the effects of five sources of uncertainty on z_orig estimates.

6.1 Measurement errors

The first source of uncertainty is measurement errors. We recalculate z_orig assuming an error of 0.4 ‰ on δD₀ (typical of what we can measure with in-situ laser instruments;  Aemisegger et al., 2012; Benetti et al., 2014) and 1 ‰ on δD_f (larger errors due to lower humidity and the increased complexity of measurements in altitude). The averaged errors on z_orig and their standard deviations are plotted as a function of EIS in Fig. 14a. Whereas errors on δD_f lead to errors on z_orig of the order of 20 m (Fig. 14a, green), errors on δD₀ lead to errors on z_orig of the order of 80 m (Fig. 14a, red). Yet in stratocumulus, no one expects the air to originate from a higher altitude than 80 m above the inversion. Therefore, δD₀ measurements would need to be more accurate than usual to be useful in stratocumulus regions, i.e., 0.1 ‰ to yield a 20 m precision on z_orig. In trade-wind cumulus regions, the precision of 0.4 ‰ is enough for z_orig to be useful.

6.2 Neglecting rain evaporation

The second source of uncertainty is associated with neglecting rain evaporation. This effect can be quantified in a model, but it is very difficult to quantify in nature because it is complicated and uncertain to measure η (Rosenfeld and Mintz, 1988), and it is even more complicated to measure or predict α_evap. Rain evaporation can have a depleting or enriching effect depending on microphysical details that are too complex to be addressed here (Graf et al., 2019). Neglecting rain evaporation leads to an error of the order of 500 m in regions of low EIS and 250 m in regions of strong EIS (Fig. 14b, brown). In regions of stratocumulus regions, rain evaporation is a significant source of error in spite of the relatively small amount of precipitation available to evaporate. This is because total evaporation of the rain efficiently enriches the SCL and easily modifies δD₀ by more than the 0.1 ‰ targeted precision explained above. However, it is possible that LMDZ overestimates this source of error in trade-wind cumulus and stratocumulus regions. LMDZ is one of the GCMs producing the strongest rain in stratocumulus regions (Zhang et al., 2013), and GCMs are known to trigger convection too often in trade-wind cumulus regions (Nuijens et al., 2015a, b).

6.3 Neglecting horizontal advection

The third source of uncertainty is associated with horizontal advection. In nature, ϕ can be estimated from meteorological analyses and β can be estimated from near-surface isotopic measurements at several locations (e.g., sounding arrays during typical field campaigns). In absence of these additional measurements, neglecting this effect leads to an error of the order of 800 m (Fig. 14b, purple). This limits the usefulness of z_orig estimates for all cloud regimes.

6.4 Daily variability in the steepness of δD profiles

The fourth source of uncertainty arises from the daily variability in α_eff (Appendix D2). Estimating α_eff requires us to measure δD_f at 500 hPa. Satellite measurements are available but are affected by random errors that are too large for our application (Worden et al., 2011, 2012; Lacour et al., 2015). Precise in situ measurements of water vapor δD in altitude are costly and difficult (Sodemann et al., 2017).

Let us assume that we have only one δD_f value that represents the seasonal average at a given location. To estimate the resulting error on z_orig, we re-estimate z_orig every day and at each location using $\stackrel{\mathrm{‾}}{{\mathit{\alpha }}_{\mathrm{eff}}}+{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{eff}}}$ and $\stackrel{\mathrm{‾}}{{\mathit{\alpha }}_{\mathrm{eff}}}-{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{eff}}}$. The error on z_orig is calculated as $\left(z_{\mathrm{orig}}(\stackrel{\mathrm{‾}}{{\mathit{\alpha }}_{\mathrm{eff}}}-{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{eff}}})-z_{\mathrm{orig}}(\stackrel{\mathrm{‾}}{{\mathit{\alpha }}_{\mathrm{eff}}}+{\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{eff}}})\right)/\mathrm{2}$. The averaged error and its standard deviation is plotted as a function of EIS in Fig. 14c (black). It is of the order of 400 m and rarely below 200 m. If we attempt to estimate α_eff as the fractionation coefficient as a function of local temperature, errors would be even more dissuasive (Fig. 14c, blue).

Therefore, estimating z_orig from daily δD₀ measurements cannot be useful unless we measure δD_f on a daily basis as well. Practically, we could imagine measuring FT properties (δD_f) at the top of a mountain while we measure δD₀ at the sea level (e.g., on islands such as Hawaii or Réunion:  Galewsky et al., 2007; Bailey et al., 2013; Guilpart et al., 2017).

6.5 Rayleigh assumption for the shape of δD profiles

Finally, as a fifth source of uncertainty comes the assumption that the δD profile follows a Rayleigh distillation line (Sect. 2.2). However, in both LMDZ (Appendix D1) and nature (Sodemann et al., 2017), δD profiles are usually intermediate between Rayleigh and mixing lines. The precision of our z_orig estimate is at a maximum in the Rayleigh distillation case.

When trying to find a numerical solution for z_orig directly from Eq. (6), a solution can be found only in 0.1 % of cases. This is because simulated δD profiles are often close to a mixing line in the lower troposphere (Appendix D1). Whatever z_orig in the lower troposphere, the δD₀ calculated from Eq. (6) is nearly constant because the δD profile is close to a mixing line (Appendix A, Fig. 4b). Whatever z_orig in the middle troposphere, the δD₀ calculated from Eq. (6) is also nearly constant because r_orig there is very small. So whatever z_orig, the δD calculated from Eq. (6) is nearly constant, and the numerical solution fails.

However, it is possible that δD profiles simulated by LMDZ are closer to mixing lines than real profiles since GCMs are known to overestimate vertical mixing through the troposphere (Risi et al., 2012b) and to mix the lower free troposphere too frequently by deep convection in trade-wind regions (Nuijens et al., 2015a, b). Therefore, the shape of δD profiles simulated by LMDZ is not a sufficient reason to reject the Rayleigh assumption. The uncertainty associated with this assumption is very difficult to quantify in LMDZ. More measurements of full δD profiles are very welcome to help quantify it.

To summarize, δD₀ measurements could potentially be useful to estimate z_orig with a useful precision, but only if we measure daily δD_f in the mid-troposphere, if the shape of δD profiles can be better documented, if we measure δD₀ at different places to quantify the effect of horizontal advection, and if we can invent innovative techniques to better quantify the effect of rain evaporation. In addition, in stratocumulus clouds, we need to measure δD₀ with an accuracy of 0.1 ‰.

Figure 14 Errors when estimating z_orig from δD₀ observations, as a function of EIS, as predicted by LMDZ. (a) Error we would make if δD₀ is measured with a 1 ‰ error (red), and error we would make if δD_f is measured with a 1 ‰ error (green). (b) Error we would make if we neglect horizontal advection effects (purple) and rain evaporation effects (brown). (c) Error if one uses α_eq as a function of local temperature to estimate α_eff (blue), error if one uses the seasonal-mean profile instead of the daily profile to estimate α_eff (black). The standard deviations among all daily errors estimated in each bin of EIS are also shown.

7 Conclusions

We propose an analytical model to predict the water vapor isotopic composition δD₀ of the sub-cloud layer (SCL) over tropical oceans. This model relies on the hypothesis that the altitude from which the air originates, z_orig, is an important factor. We build on B15, who extended the Merlivat and Jouzel (1979) closure equation to make the explicit link between δD₀ and mixing processes. We further extend their equation: we assume a shape for the δD vertical profiles as a function of q, and we account for horizontal advection and rain evaporation effects.

The resulting equation highlights the fact that δD₀ is not sensitive to the intensity of mixing processes. Therefore, it is unlikely that water vapor isotopic measurements could help estimate the entrainment velocity that many studies have striven to estimate (Bretherton et al., 1995). In contrast, δD₀ is sensitive to the altitude from which the air originates. Based on a simulation with LMDZ and observations during the STRASSE cruise, we show that z_orig is an important factor explaining the seasonal spatial and daily variations in δD₀, especially in subsidence regions. In turn, could δD₀ measurements, combined with vertical profiles of humidity, temperature, and δD, help estimate z_orig and thus discriminate between different mixing processes? For such isotope-based estimates of z_orig to be useful, we would need a precision of a few hundreds of meters in deep convective regions and smaller than 20 m in stratocumulus regions. To reach this target, we would need daily measurements of δD in the mid-troposphere and very accurate measurements of δD₀, which are currently difficult to obtain. We would also need information on the horizontal distribution of δD to account for horizontal advection effects, and full δD profiles to quantify the uncertainty associated with the assumed shape for δD profiles. Finally, rain evaporation is an issue in all regimes, even for stratocumulus clouds. Innovative techniques would need to be developed to quantify this effect from observations.

This study is preliminary in many respects. First, it would be safe to check using water tagging experiments in LMDZ that z_orig estimates really represent the altitude from which the air originates and are not biased by our simplifying assumptions. Second, the coarse vertical resolution of LMDZ and the simplicity of mixing parameterizations (e.g., cloud top entrainment is not represented) are a limitation of this study. Ideally, the relationship between δD₀, z_orig, and the type of mixing processes should be investigated in isotope-enabled large-eddy simulations (LESs) (Blossey et al., 2010; Moore et al., 2014). Artificial tracers and structure detection methods (Park et al., 2016; Brient et al., 2019), combined with conditional sampling methods (Couvreux et al., 2010), could help detect the different kinds of mixing structures, estimate their contributions to vertical transport, and describe their isotopic signature. This would allow us to confirm, or disprove, many of the hypotheses and conclusions in this paper. Finally, if the sensitivity of δD₀ to the type of mixing processes is confirmed, paired isotopic simulations of single-column model (SCM) versions of general circulation models (GCMs) and LES, forced by the same forcing, could be very useful to help evaluate and improve the representation of mixing and entrainment processes in GCMs, as is routinely the case for non-isotopic variables (Randall et al., 2003; Hourdin et al., 2013; Zhang et al., 2013).

Appendix A: Closure if the tropospheric profile follows a mixing line

For simplicity, we neglect here horizontal advection and rain evaporation effects, but results would be similar otherwise. If we assume that R_orig is uniquely related to q_orig through a mixing line between the SCL air and a dry end member (q_f, R_f),

$$\begin{array}{}\text{(A1)}& q_{\mathrm{orig}}=a\cdot q_{\mathrm{0}}+(\mathrm{1}-a)\cdot q_{\mathrm{f}}\end{array}$$

and

$$\begin{array}{}\text{(A2)}& R_{\mathrm{orig}}=a\cdot q_{\mathrm{0}}\cdot R_{\mathrm{0}}+(\mathrm{1}-a)\cdot q_{\mathrm{f}}\cdot R_{\mathrm{f}}.\end{array}$$

Reorganizing Eq. (A1), we get $a=\frac{r-p}{\mathrm{1}-p}$ with $p=q_{\mathrm{f}}/q_{\mathrm{0}}$. Since $q_{\mathrm{f}}\le q_{\mathrm{orig}}\le q_{\mathrm{0}}$, p≤r_orig. Injecting Eq. (A2) into Eq. (6), we get

$$\begin{array}{}\text{(A3)}& R_{\mathrm{0}}={\displaystyle \frac{R_{\mathrm{oce}}/{\mathit{\alpha }}_{\mathrm{eq}}+p/(\mathrm{1}-p)\cdot R_{\mathrm{f}}\cdot {\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})}{h_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{K}}\cdot (\mathrm{1}-h_{\mathrm{0}})/(\mathrm{1}-p)}}.\end{array}$$

As a consistency check, in the limit case where the end member is totally dry (p=0), we find the MJ79 equation, i.e., Eq. (10).

It is intriguing to realize that r_orig has disappeared from Eq. (A3). This can be understood physically: if the vertical profile follows a mixing line, it does not matter from which altitude the air comes: ultimately, what matters is how much dry air has been mixed directly or indirectly into the SCL (Fig. 4b). Therefore, if R_orig follows a mixing line, we lose the sensitivity to z_orig.

Appendix B: Diagnostics for rain evaporation in LMDZ

Rain evaporation can be accounted for in Eq. (8) if we can quantify η, the ratio of water vapor originating from rain evaporation to that originating from surface evaporation, and α_evap, the ratio of isotopic ratio in the rain evaporation flux to R₀.

B1 Equations

In LMDZ, two parameterization schemes can produce rain evaporation: the convective scheme and the large-scale condensation scheme. Their respective precipitation evaporation tendencies, ${\left(\mathrm{d}q/\mathrm{d}t\right)}_{\mathrm{evap},\mathrm{conv}}$ and ${\left(\mathrm{d}q/\mathrm{d}t\right)}_{\mathrm{evap},\mathrm{lsc}}$, are given in ${\mathrm{kg}}_{\mathrm{water}}\cdot {{\mathrm{kg}}_{\mathrm{air}}}^{-\mathrm{1}}\cdot {\mathrm{s}}^{-\mathrm{1}}$ and are used to calculate F_evap in ${\mathrm{kg}}_{\mathrm{water}}\cdot {\mathrm{m}}^{-\mathrm{2}}\cdot {\mathrm{s}}^{-\mathrm{1}}$:

$$F_{\mathrm{evap}}=\sum _{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left({\left(\mathrm{d}q/\mathrm{d}t\right)}_{\mathrm{evap},\mathrm{conv}}+{\left(\mathrm{d}q/\mathrm{d}t\right)}_{\mathrm{evap},\mathrm{lsc}}\right)\cdot {\displaystyle \frac{\mathrm{\Delta }{P}_k}{g}},$$

where k_LCL is the last layer below the LCL, ΔP_k is the depth of layer k in pressure coordinate, and g is gravity.

The isotopic equivalent of this flux, F_evap,iso, is used to calculate $R_{\mathrm{evap}}=F_{\mathrm{evap},\mathrm{iso}}/F_{\mathrm{evap}}$.

Only grid boxes and days where F_evap>0.05 mm d⁻¹ are considered to calculate α_evap. This represents 94.0 % of all tropical oceanic grid boxes.

B2 Results

Consistent with the larger amount of precipitation available for evaporation, η is at a maximum in regions of deep convection, reaching 30 % around the Maritime Continent (Fig. C1a). It is minimal over the dry descending regions, reaching 5 % off the coasts of Mauritania, Peru, and Namibia. The rain evaporation is more depleted than the SCL in regions of strong deep convection, by as much as 70 ‰ around the Maritime Continent (Fig. C1b). When the fraction of raindrops that evaporate is small, as is the case in such moist regions, isotopic fractionation favors evaporation of the lighter isotopologues. In these regions, rain evaporation has a depleting effect on the SCL, consistent with Worden et al. (2007). In contrast, in other regions, rain evaporation has an enriching effect on the SCL, up to 70 ‰ in dry regions. This is because in dry regions, rain evaporates almost totally, so that the evaporation flux has almost the same composition as the initial rain, which is more enriched than the water vapor.

Appendix C: Diagnostics for horizontal advection in LMDZ

We can account for horizontal advection in Eq. (8) if we can quantify parameters $\mathit{\varphi }=\frac{F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}}{E}$, the ratio of water vapor coming from horizontal advection to that coming from surface evaporation, and $\mathit{\beta }=R_{\mathrm{adv}}/R_{\mathrm{0}}$, the ratio of isotopic ratios of horizontal advection to those of the SCL.

C1 Equations

Let us assume that the box representing the SCL has a zonal extent Δy and a meridional extent Δx and is composed of k_LCL layers of vertical extent Δz_k. The quantity F_adv⋅q_adv represents the mass flux of water entering the grid box by horizontal advection per surface area, expressed in ${\mathrm{kg}}_{\mathrm{water}}\cdot {\mathrm{s}}^{-\mathrm{1}}\cdot {\mathrm{m}}^{-\mathrm{2}}$. Assuming an upstream advection scheme, it can be expressed as

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}=\\ & {\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left({\mathit{\rho }}_k\cdot \left|u_k\right|\cdot q_{uk}\cdot \mathrm{\Delta }y\cdot \mathrm{\Delta }{z}_k+{\mathit{\rho }}_k\cdot \left|v_v\right|\cdot q_{vk}\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }{z}_k\right)}{\mathrm{\Delta }x\cdot \mathrm{\Delta }y}}\end{array},\end{array}$$

where u_k and v_k are the zonal and meridional wind components at layer k, ρ_k is the volumic mass of air at layer k, and q_uk and q_vk are the humidities of the incoming air from zonal and meridional advection at layer k. When u_k>0, q_uk is the humidity in the grid box to the west. When u_k<0, q_uk is the humidity in the grid box to the east. When v_k>0, q_vk is the humidity in the grid box to the south. When v_k<0, q_vk is the humidity in the grid box to the north.

Applying the hydrostatic equation at each layer ($\mathrm{\Delta }{P}_k={\mathit{\rho }}_k\cdot g\cdot \mathrm{\Delta }{z}_k$, where g is gravity and ΔP_k is the vertical extent of the layer k in pressure coordinate), we get

$$F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}=\sum _{k=\mathrm{1}}^{k_{\mathrm{LCL}}}{\displaystyle \frac{\mathrm{\Delta }{P}_k}{g}}\cdot \left({\displaystyle \frac{\left|u_k\right|\cdot q_{uk}}{\mathrm{\Delta }x}}+{\displaystyle \frac{\left|v_k\right|\cdot q_{vk}}{\mathrm{\Delta }y}}\right).$$

The quantity F_adv represents the incoming air mass flux by horizontal advection, and q_adv represents the humidity of the incoming air. We can thus write them as

$$F_{\mathrm{adv}}=\sum _{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left({\displaystyle \frac{\left|u_k\right|}{\mathrm{\Delta }x}}+{\displaystyle \frac{\left|v_k\right|}{\mathrm{\Delta }y}}\right)\cdot {\displaystyle \frac{\mathrm{\Delta }{P}_k}{g}}$$

and

$$q_{\mathrm{adv}}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left(\frac{\left|u_k\right|}{\mathrm{\Delta }x}\cdot q_{uk}+\frac{\left|v_k\right|}{\mathrm{\Delta }y}\cdot q_{vk}\right)\cdot \frac{\mathrm{\Delta }{P}_k}{g}}{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left(\frac{\left|u_k\right|}{\mathrm{\Delta }x}+\frac{\left|v_k\right|}{\mathrm{\Delta }y}\right)\cdot \frac{\mathrm{\Delta }{P}_k}{g}}}.$$

The same budget as in Eq. (C1) can be written for water isotopes:

$$\begin{array}{rl}& F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}\cdot R_{\mathrm{adv}}=\\ & {\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left({\mathit{\rho }}_k\cdot \left|u_k\right|\cdot q_{uk}\cdot R_{uk}\cdot \mathrm{\Delta }y\cdot \mathrm{\Delta }{z}_k+{\mathit{\rho }}_k\cdot \left|v_v\right|\cdot q_{vk}\cdot R_{vk}\cdot \mathrm{\Delta }x\cdot \mathrm{\Delta }{z}_k\right)}{\mathrm{\Delta }x\cdot \mathrm{\Delta }y}}\end{array},$$

where R_adv represents the isotopic ratio of the incoming water vapor:

$$R_{\mathrm{adv}}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left(\frac{\left|u_k\right|}{\mathrm{\Delta }x}\cdot q_{uk}\cdot R_{uk}+\frac{\left|v_k\right|}{\mathrm{\Delta }y}\cdot q_{vk}\cdot R_{vk}\right)\cdot \frac{\mathrm{\Delta }{P}_k}{g}}{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}\left(\frac{\left|u_k\right|}{\mathrm{\Delta }x}\cdot q_{uk}+\frac{\left|v_k\right|}{\mathrm{\Delta }y}\cdot q_{vk}\right)\cdot \frac{\mathrm{\Delta }{P}_k}{g}}}.$$

Note that the upstream advection scheme assumed here overestimates the effect of advection compared to the Van Leer (1977) advection scheme used in LMDZ. We thus estimate an upper bound for the advection effect here.

In practice, rather than calculating $\mathit{\beta }=R_{\mathrm{adv}}/R_{\mathrm{0}}$ , we calculate $\mathit{\beta }=R_{\mathrm{adv}}/R_{\mathrm{SCL}}$, where R_SCL is the isotopic ratio on average through the SCL:

$$R_{\mathrm{SCL}}={\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}{q}_k\cdot R_k\frac{\mathrm{\Delta }{P}_k}{g}}{{\sum }_{k=\mathrm{1}}^{k_{\mathrm{LCL}}}{q}_k\frac{\mathrm{\Delta }{P}_k}{g}}}.$$

This prevents the advected water vapor to be systematically more depleted when the mixed-layer hypothesis is not exactly verified.

C2 Results

Parameter ϕ is at a maximum where winds are maximum, such as near the extra-tropics or in the North Atlantic (Fig. C1a). Horizontal advection has an enriching effect in deep convective regions (probably because water vapor comes from nearby drier regions that have been less depleted by deep convection) and a depleting effect near the coasts (probably because of winds bringing vapor from the nearby land that is depleted by the continental effect) (Fig. C1b).

Note that in this formulation, parameters ϕ and β are resolution-dependent. For example, in a finer resolution, ϕ would be larger and β would be closer to 1, but $F_{\mathrm{adv}}\cdot q_{\mathrm{adv}}\cdot R_{\mathrm{adv}}$ and thus the contribution of horizontal advection in Eq. (8) would remain the same.

Figure C1 (a) Ratio η of water in the SCL coming from rain evaporation to that coming from surface evaporation. (b) Effective fractionation coefficient α_evap between the SCL water vapor and the rain evaporation. (c) Ratio ϕ of water in the SCL coming from horizontal advection to that coming from surface evaporation. (d) Effective fractionation coefficient β between the SCL water vapor and the water vapor coming from horizontal advection.

Appendix D: LMDZ free-tropospheric profiles

The goal of this appendix is to document the spatiotemporal variability in the shape (Sect. D1) and steepness (Sect. D2) of simulated free-tropospheric δD profiles. Note that a detailed interpretation of these profiles is beyond the scope of this paper. This paper aims at understanding δD₀, which is the first step towards understanding full tropospheric profiles. In turn, understanding full tropospheric profiles in future studies will help refine our model for δD₀.

Figure D1 (a) Vertical profiles of water vapor δD simulated by LMDZ (black), calculated if assuming a Rayleigh-type curve with α_eff estimated from (q_f,δD_f) at 500 hPa (green), calculated if assuming a Rayleigh-type curve with α_eq(T) (dashed green), and calculated if assuming a mixing curve between the first layer and (q_f,δD_f) at 500 hPa (blue), on average over all tropical oceanic locations and days. (b) Winter-mean map of parameter $f=\frac{\mathit{\delta }{D}_{\mathrm{LMDZ}}-\mathit{\delta }{D}_{\mathrm{Rayleigh}}}{\mathit{\delta }{D}_{mix}-\mathit{\delta }{D}_{\mathrm{Rayleigh}}}$ at 4000 m, i.e., slightly below 500 hPa where q_f and δD_f are taken. Parameter f describes how close the simulated δD (δD_LMDZ) is to the Rayleigh (δD_Rayleigh) and mixing lines (δD_mix): f=0 in case of a Raleigh line, f=1 in case of a mixing line, and f>1 if δD is more enriched than a mixing line. (c) Standard deviation of parameter f among all days in winter of all years, at 4000 m. (d) Same as (b) but at 1000 m, i.e., slightly above the LCL. (e) Same as (c) but at 1000 m. To avoid numerical problems, only days and locations where $|\mathit{\delta }{D}_{\mathrm{mix}}-\mathit{\delta }{D}_{\mathrm{Rayleigh}}|>\mathrm{5}$ ‰ are used in the calculations.

D1 Shape of tropospheric profiles

First, we test whether the δD vertical profiles simulated by LMDZ follow a Rayleigh or mixing line as a function of q. For the Rayleigh curve, α_eff is estimated as explained in Sect. 3.3. For the mixing line (Appendix A), the end member (q_f, R_f) is also taken at 500 hPa.

The tropical-mean vertical δD profiles simulated by LMDZ are bounded by Rayleigh and mixing lines (Fig. D1a). To better document the spatial variability in the shape of δD profiles, we plot parameter $f=\frac{\mathit{\delta }{D}_{\mathrm{LMDZ}}-\mathit{\delta }{D}_{\mathrm{Rayleigh}}}{\mathit{\delta }{D}_{\mathrm{mix}}-\mathit{\delta }{D}_{\mathrm{Rayleigh}}}$, describing how close the simulated δD (δD_LMDZ) is to the Rayleigh (δD_Rayleigh) and mixing (δD_mix) lines. We have f=0 in the case of a Raleigh line, f=1 in the case of a mixing line, and f>1 if δD is more enriched than a mixing line. In the lower troposphere, δD_LMDZ is close to a mixing line (and sometimes even more enriched) in deep convective regions (Indian Ocean, South Pacific Convergence Zone, Atlantic ITCZ), probably because deep convection efficiently mixes the lower troposphere. Elsewhere, δD_LMDZ is intermediate between the two lines (Fig. D1e). In the middle troposphere, δD_LMDZ is relatively closer to Rayleigh everywhere (Fig. D1b).

The daily variability of f is large everywhere and at all levels (Fig. D1c, e), with standard deviation of 0.23 and 0.44 on tropical average at 1000 and 4000 m, respectively. A large daily variability in the shape of profiles is also observed in nature (Sodemann et al., 2017).

D2 Steepness of tropospheric profiles

The steepness of the δD gradient from the surface to the middle troposphere is described by the parameter α_eff. It is at a maximum in regions of deep convection, for example around the Maritime Continent (Fig. D2a). This is consistent with the maximum depletion simulated in deep convective regions in the mid-troposphere simulated by models (Bony et al., 2008), leading to steeper δD profiles. The pattern of α_eff may also reflect horizontal advection effects, where strong isotopic gradients align with winds (e.g., from the eastern to the western Pacific;  Dee et al., 2018).

Values of α_eff are of the same order of magnitude as real fractionation factors, but the spatial variations do not reflect those predicted if using a fractionation coefficient α_eq as a function of temperature T (Fig. D2b).

The daily standard deviation of α_eff (${\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathrm{eff}}}$) for a given season ranges from 5 ‰ in the central Atlantic to 40‰ near the Maritime Continent (Fig. D2c). On average over all seasons and locations, daily α_eff−1 at a given location varies within ±25 % of its seasonal-mean mean value.

Figure D2 (a) α_eff−1, where α_eff is the effective fractionation coefficient, expressed in per mill. (b) α_eq(T)−1 expressed in per mill. All daily values are averaged over all days in winters of all years. (c) Standard deviation of α_eff among all days in winter of all years, expressed in per mill.