Beyond Craig and Gordon: A model of water vapor isotopologues in the marine boundary layer

We develop a one-dimensional (1D) steady state isotope marine boundary layer (MBL) model that includes meteorologically important features missing in Craig and Gordon type models, namely height-dependent diffusion/mixing, lifting to deliver air to the free troposphere, and convergence of subsiding air. Kinetic isotopic fractionation results from this 10 height-dependent diffusion that starts as pure molecular diffusion at the air-water interface and increases with height due to turbulent eddies. Convergence causes mixing of dry, isotopically depleted air with ambient air. Model results fill a quadrilateral in δD-δO space, of which three boundaries are respectively defined by 1) vapor in equilibrium with various sea surface temperatures (SSTs); 2) mixing of vapor in equilibrium with seawater and vapor in subsiding air; and 3) vapor that has experienced maximum possible kinetic fractionation. Model processes also cause variations in d-excess of MBL 15 vapor. In particular, mixing of relatively high d-excess descending/converging air into the MBL increases d-excess, even without kinetic isotope fractionation. The model is tested by comparison with seven datasets of marine vapor isotopic ratios, with excellent correspondence. About 95% of observational data fall within the quadrilateral predicted by the model. The distribution of observations also highlights the significant influence of vapor from nearby converging descending air on isotopic variations within the MBL. At least three factors may explain the ~5% of observations that fall slightly outside of 20 the predicted regions in δD-δO and d-excess-δO space: 1) variations in seawater isotopic ratios, 2) variations in isotopic composition of subsiding air, and 3) influence of sea spray.


Introduction
Stable isotopic ratios of water have been widely used to study the hydrologic cycle of the atmosphere.They have proven to be a powerful tool for understanding modern atmospheric processes (e.g., Dansgaard, 1964, Lawrence et al., 2004; Worden 25 et al., 2007;Uemura et al., 2008;Kurita, 2011;Kopec et al., 2017).In addition, they have been extremely useful for inferring paleoclimate conditions and making climate reconstructions from glacier ice, tree rings, lake sediments, speleothems, and paleosols (e.g., Dansgaard et al., 1989;Wang et al., 2001;Huang et al., 2002;Andersen et al., 2004;Jouzel et al., 2007;Feng et al., 2007;Sheldon and Tabor, 2009;Masson-Delmotte et al., 2015).Particularly relevant to this study is the adaptation of the C-G model for the marine boundary layer.An influential study by Merlivat and Jouzel (1979) linked the magnitude of kinetic isotopic fractionation primarily within the laminar layer above the water-air interface, which is required input for the C-G model, to aerodynamic conditions, i.e., wind speed and surface ocean roughness.However, the model still required the input of the free atmosphere humidity and isotopic ratios.
Recognizing the difficulty of knowing these free atmosphere variables, Merlivat and Jouzel (1979) made an assumption, 5 known later as the "closure assumption", that the isotopic ratios of vapor mass in the free atmosphere were equal to the isotopic ratios of the vapor fluxes from the sea surface.This assumption enabled them to complete a new multi-slab model (the M-J model), used by numerous investigators to calculate isotopic fluxes from the sea surface over a range of maritime conditions, and to explore relationships between isotopic compositions of evaporative flux and boundary layer meteorological conditions such as sea surface temperature and relative humidity (e.g., Johnsen et al., 1989;Petit et al. 1991).10 The closure assumption also allowed the modeled flux to be used as the starting isotopic composition of an air mass, which evolves during transport and rainout or a Rayleigh process (e.g., Johnsen et al., 1989;Petit et al. 1991).The closure assumption has been determined to be generally invalid at local scales (Jouzel and Koster, 1996), however, it has continued to be used (e.g., Benetti et al., 2014) simply for lack of a better assumption.
Abandoning this closure assumption requires a fundamental rethinking of the MBL model structure.In addition, there are 15 ramifications for other model assumptions.As a consequence, we consider three requirements for developing a physically consistent MBL model free of the invalid closure assumption.
First, vertical advection is necessary at the inception of a Rayleigh process in order to lift MBL air into the free atmosphere.
When an air mass is lifted into the free troposphere, the vapor isotopic ratio of the air for the first condensation is equal to the isotopic ratio of vapor within the air mass, not the ratio of isotopologue diffusive fluxes into the air mass.Contrary to the 20 closure assumption, these two are not generally equal.Therefore, an MBL model should calculate not only the isotopic ratio of vapor flux, but also that of vapor within the MBL, and the latter is the quantity that should be used for the initial vapor isotopic composition in any subsequent Rayleigh process.
Second, with incorporation of vertical advection, mass balance requires 1) horizontal convergence of air within the MBL to replenish the lifted air in the evaporation column, and 2) subsidence of air outside the model region to sustain the local 25 horizontal convergence.Such a circulation on various scales was discussed by Craig and Gordon (1965) to explain why vapor in the MBL was not in isotopic equilibrium with ocean water.In this contribution, we attempt to quantify how horizontal convergence of air from non-local regions of subsidence affects the isotopic properties of the local MBL.Because the converging air is unlikely to have the same isotopic composition as the local MBL air, convergence turns out to affect MBL vapor isotopes quite significantly, as discussed later in this paper . 30 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
Third, incorporating convergence means that the assumption of constant flux in the C-G model is no longer valid, even under steady state conditions.This necessitates different equations of mass conservation.
In addressing these three required changes, we incorporate a fourth major change to the model structure, which is to convert previous multi-slab models (referred to hereafter as C-G type models) to a true one-dimensional model.In doing so, we remove the laminar slab near the sea surface and describe the kinetic fractionation as a height dependent process that is 5 controlled by increasing turbulent transport coefficient with height (see below).An additional benefit is the ability to obtain isotopic ratios of air and vapor flux at any given height within the MBL.
Because advective and diffusive fluxes vary with height in a 1D model, careful choice of parameterization of the diffusive fluxes can minimize additional complexity.The eddy diffusion coefficient (coefficient of turbulent transport in eddies) increases continuously with height (Merlivat and Coantic, 1975) from zero at the air-water interface, where transport of 10 vapor is thus effected solely by isotopically fractionating molecular diffusion, to greater values at height where vapor transport is by turbulent eddies.Such a height-dependent change of diffusion coefficient is adopted in our model.In addition, our model is relieved of the need 1) to empirically choose the value of a parameterized kinetic fractionation factor (Δε in Craig and Gordon, 1965;k in Merlivat and Jouzel, 1979) (this choice may sometimes be difficult and values reported in literature may not apply to specific conditions under investigation (e.g., Xiao et al., 2017)), and 2) to specify a specific 15 laminar layer thickness.Instead, we allow the diffusion coefficient to represent pure molecular diffusion at the interface (which differs for different isotopologues and thus leads to fractionation), and to increase linearly with height to several orders of magnitude greater than molecular diffusivities.
Abandoning some of the assumptions of earlier models, such as constant flux of vapor isotopologues, flux ratios equal to concentration ratios (the invalid closure assumption), and presence of a discrete laminar layer, permits a significantly more 20 realistic and elucidative approach to understanding processes in the MBL and allows more meteorological profiles of variables (such as humidity and isotopic ratios of the vapor in the MBL) to be calculated rather than specified.The tradeoff, obviously, is in sacrificing the simplicity of the classical model.The model reported here attempts to balance that tradeoff: it is considerably less complex than isotope-enabled general circulation models (GCMs) and should be accessible to investigators without substantial experience with complex models, yet it allows exploration of physical controls of vapor 25 within the MBL and in the initial Rayleigh process above the MBL.understanding the systematics linking isotopic observations of precipitation to the meteorological conditions of the vapor source, of the precipitation site, and along the moisture paths between the two.We use the new MBL model presented here to examine the vapor source part of the isotope systematics.Since the model produces vapor concentrations and isotopic ratios, it can be tested and validated by MBL isotopic measurements, which, thanks to new spectral vapor isotopic measurement technology, have become increasingly available.There are still additional potential benefits.For example, 10 such a model might provide a new way to estimate evaporation rate, one of the holy grails of weather and climate models.
In the following sections, we first describe the formulation and solution of the model and the marine boundary layer observations to be used to validate the model.Then we discuss the model results and their comparison with the observations, as a basis for addressing the systematics of vapor source conditions and atmospheric isotopes.Although the limitations of the model will be discussed in more detail in section 6.2, we briefly mention here that this model applies to the part of the 15 marine boundary layer where vertical velocity is positive (upward), there is no net horizontal advection, and the model does not include vapor liquid exchange within the air column.

The Isotope Marine Boundary Layer Model
The model we describe here has been developed to study the effect of marine boundary layer processes, such as evaporation of water, mixing and uplift of air, on concentrations and fluxes of isotopologues of the MBL.Three isotopologues, H 2 16 O, 20 H 2 18 O and HDO are modeled and presented here, but more can be added easily.We refer to this model as the Isotope Marine Boundary Layer (IMBL) model.
Figure 1 is a cartoon of the IMBL model showing the three layers that comprise the model column itself, and the input of external air.Layer 1, the lowest layer, extending from the surface at z=0 to height z=h 1 , is a quasi-von Kármán layer in which vapor is transported upward from the sea surface by mixing that increases in intensity with height.Layer 2 (h 1 <z<h 2 ), 25 the middle layer, is subject to strong vertical mixing, to the convergence of air that has elsewhere descended from the free atmosphere and converged horizontally into the modeled column, and to vertical advection caused by the convergence.In Layer 3 (h 2 <z<h 3 ), the top layer, there is no convergence, so the air ascends at a fixed rate, while the vertical mixing rate decreases in intensity with height.Brien (1970).The coefficient K(z) equals the molecular diffusion coefficient K m at the surface and increases linearly with height to a maximum value K max at z=h 1 .It remains fixed at K max through the middle layer, then decreases linearly in the top layer above z=h 2 to a small value K t at z=h 3 .The vertical velocity w(z) is zero in 5 Layer 1, increases linearly with height through the middle layer, in which convergence occurs at a fixed rate, and remains constant at value w a in Layer 3. Consistent with their constant values of w, Layers 1 and 3 do not have convergence.
The following subsections, 2.1-2.3,describe the individual physical and mathematical features of the model.Table 1 contains a list of variables and parameters found in these subsections and elsewhere.

Mixing Process 10
The central matter for this subsection is the specification of the height-dependent eddy diffusion coefficient,  !(z), which appears in Fick's Law for diffusive flux, where  ! is the vertical flux of the i'th isotopologue (isotopologue-mass area -1 time -1 ), C i is the concentration ratio of the i'th isotopologue (isotopologue-mass dry-air-mass -1 ),  is the density of dry air (mass-of-dry-air volume -1 ), and  is the vertical 15 coordinate (increasing upwards from z=0 at the surface).The i-subscripts of F, K and C are reminders that they all depend on the specific isotopologue under consideration, but for simplicity we drop them hereafter.We note that C i has the same units as the commonly used mixing ratio r.The difference is that r is total water vapor mass per unit dry air mass, while C i is the mass of the i'th isotopologue (e.g., H 2 18 O) per unit dry air mass.In this paper, we will use the term concentration ratio for C and mixing ratio for r.With (1), we assume that Fick's Law can be used to represent vertical mixing by the combined effects 20 of mechanically-driven turbulence, buoyancy-driven convection, and molecular diffusion.
In adopting Fick's Law, here, we have made the tacit assumption that alternative mixing models are less appropriate for our purposes.While higher-order closure schemes (e.g., Burk, 1977), structured turbulence models (e.g., Kirwan, 1968), and the telegraph equation (e.g., Goldstein, 1951) have some advantages over Fick's Law, their added complexity would not be justified at this juncture, and we postpone their consideration until future investigations warrant.25 . (2) For  given by Eq. ( 1), and for the case of  with negligible dependence on  or , Eq. (2) becomes

𝐾(𝑧)
!" !" . (3) Returning to the central matter, the specification of (z), we reject the assumption of constant , the simplest and most frequently used assumption, because it is particularly unrealistic near boundaries (i.e., water-air interface in this work), 5 where the inhibitive effect of the interface on mixing of air increases with proximity to the boundary.The next most frequently used assumption is that  is a linear function of , although a few others have been proposed (Merlivat and Coantic, 1975).The use of linear functions of z to represent   has a long history in turbulence studies, including the turbulent transport of momentum as well as both buoyantly active and passive scalar fluid properties.The well-known work of von Kármán (1930) andPrandtl (1932) successfully applied the simple form   =  • , where b is a constant, to derive 10 the equation of the logarithmic layer, where   =  • ln  + , with u being the wind speed, and s and m being constants.
An obvious limitation of the widely cited von Kármán/Prandtl formulation occurs when z is very small, near the singularity at z=0.The most common way to circumvent this problem has been to introduce a discrete "laminar boundary layer" (LBL), a very thin but finite layer with constant diffusion by molecular motion and with weak turbulent influence.The incremental cost of this approach is the necessity of specifying one additional parameter, δ, the thickness of the LBL.15 Another way to overcome the problem for small z is to use the more general form: where  ! is the molecular diffusion coefficient for vapor in air and  •  is the contribution of turbulent eddies to the diffusion coefficient.Note that  !varies among isotopologues, but b does not.An equivalent general linear form was applied to boundary layer mixing above the LBL by Montgomery (1940), and within the LBL by Sverdrup (1946;1951).20 One advantage of the form of Eq. ( 4) for the parameterization is the gain of one degree of freedom through the use of the known quantity  !instead of the unknown parameter δ, the thickness of the laminar layer.The latter can be replaced by a diagnostic laminar layer thickness, z * , the height where molecular and turbulent diffusion coefficients are equal.In other words, below z * vertical diffusion is dominated by molecular processes, and above it turbulence and convection prevail.(5) A linear approach, mathematically equivalent to Eq. ( 4) (Sheppard, 1958), used bulk aerodynamic theory to modify Sverdrup's (1946Sverdrup's ( , 1951) ) work.The result was another linear function of z containing the friction velocity u * and von Kármán's constant κ instead of the coefficients  ! , thus connecting Sheppard's model to familiar parameters of fluid mechanics.5 Merlivat and Coantic (1975) tested and compared various linear and nonlinear alternatives to Eq. ( 4).In contrasting alternative boundary layer models for use in isotope studies, they concluded that their laboratory experiments did not support Sheppard's linear theory.However, Eq. ( 4), which is mathematically equivalent to Sheppard's (1958) approach, was used successfully to model atmospheric vapor in field experiments above Arctic lakes (Feng et al, 2016).Hence, we proceed with Eq. (4). 10

Convergence and Vertical Advection
Moist air undergoing the Rayleigh distillation process in the free atmosphere is generally conceived to have originated in the PBL and been lifted (i.e., vertically advected) into the free atmosphere.For mass to be conserved, such uplift must be accompanied by convergence within the PBL.Nevertheless, C-G type models ignore convergence within the boundary layer (e.g., Craig and Gordon, 1965;Merlivat and Jouzel, 1979).The incorporation of this apparent contradiction into a model 15 might be justified by arguing that the effect of boundary layer convergence on isotopic processes is negligible, or if the only concern is the isotopic evolution of the liquid where the vapor originates.In the IMBL model presented here, however, we choose to preserve consistency by including both convergence and uplift, and to use model results to diagnose the importance of the convergence effect rather than neglecting it a priori.As we later show, convergence has a large influence on the isotopic composition of the air exiting the MBL upward into the free atmosphere.20 Steady-state conservation of mass for dry air, using dynamic variables and neglecting diffusion, can be written in the form where  (here considered independent of ) is the horizontal dynamic convergence (dry-air-mass volume -1 time -1 ), and () is the dynamic vertical velocity (dry-air-mass area -1 time -1 ).The kinematic (conventional) velocity (length time -1 ) is the dynamic velocity divided by the air density, ρ.We will use Eq. ( 6) to derive the governing equation for the middle layer 25 in section 2.3.2.Ignoring (for now) the effect of diffusion, conservation of mass for isotopologues affected only by kinematics can now be expressed as where  ! is the concentration ratio of the isotopologue of the MBL air converging into the area being modeled.The first and second terms on the right are the direct effect of convergence (replacement of air of concentration ratio C by converging air 5 of concentration ratio C C ) and the effect of vertical advection, respectively.Note that Eq. ( 7) is also consistent with the assumed absence of non-divergent horizontal advection.
The converging air, with concentration ratio  !, is a mixture of two air types, with fractional presence by mass β and (1-β), respectively: 1) air from aloft, originally with concentration ratio  !, that has been recently integrated into the MBL by sinking or mixing, and 2) air that has been in the MBL for considerable time and has become essentially identical in 10 properties to the modeled air with concentration ratio C. Thus, C C = βC E + (1-β)C and Eq. ( 7) can be written in terms of  !as: . (8)

Governing Equations
To find the general form of the steady-state equation of conservation of mass for each vapor isotopologue, we combine the 15 diffusive (Eq. 3) and kinematic (Eq.8) effects and set (Since dynamic variables are used here, this result does not depend on the commonly invoked isopycnal approximation.) Eq. ( 9) is the general form of the basic governing equation that we solve in layers in which () and   change.This governing equation is implemented three times, once for each isotopologue, with  differing among isotopologues.20 Equivalently, it may be viewed as a single vector equation of length 3, with each component describing mass conservation for one isotopologue.The method devised here to solve Eq. ( 9), described in 3.2, uses the latter strategy.
We now proceed to adapt Eq. ( 9) to the atmospheric conditions specific to layers 1-3 --the low, middle and high layers --of the MBL.

Low layer equation
In the low layer, as described at the beginning of this section and as illustrated in Figure 1, there is no convergence or uplift.
Hence  = 0 and   = 0.As specified by Eq. ( 4),  increases linearly with height, from the small molecular value  ! at the surface ( = 0) to the larger mixing rate K max at z=h 1 , the top of the low (von Kármán) layer, where In the low layer, Eq. ( 9) thus simplifies to:

Middle layer equation
In the middle layer, where ℎ !≤  ≤ ℎ !, K(z) is constant and equal to  !"# , and the convergence rate  is also constant.
Defining w a as the upward velocity at the top of the middle layer, h 2 , Eq. ( 6) implies that w(z) increases linearly from w=0 at 10 z=h 1 to w=w a at z=h 2 .i.e., Eq. ( 9), after substituting Eqs. ( 12)-( 13), simplifies to: Within the middle layer, vertical mixing (the first term in Eq. 14) is controlled by the constant eddy diffusion coefficient,  !"# , which is the maximum value of K.The second term in Eq. ( 14) describes the direct effect of convergence of external air from aloft, originally of concentration ratio  !, into the profile.Vertical advection (the third term) occurs at a rate depending on the linearly increasing velocity and the gradient of .

High layer equation 20
The upper layer of the MBL, just below the very top, is often capped by a stable inversion in which diffusion plays a minimal role.Uplift, however, continues upward unabated through the inversion into the free atmosphere, where further evolution of the air mass is beyond the scope of the IMBL model.In the upper layer of the MBL, we assume that K(z) Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.decreases linearly from  !"# at z=ℎ ! to  ! at the top of the MBL (z=ℎ ! ), and that there is no further convergence, so   here equals w a .Eq. ( 9) thus becomes: 3 Solution Methods

Analytic Solutions 5
All three governing equations, (Eqs. 11,14 and 15), are second order linear ordinary differential equations with non-constant coefficients.Eqs. ( 11) and ( 15) are homogeneous, while Eq. ( 14) is inhomogeneous by virtue of  ! .Each equation has an analytic solution with two constants of integration, totaling six constants requiring six boundary conditions (BC's).The six BC's are: 0) is in equilibrium with the surface water.
(BC1) 10 C(z) and KdC/dz are continuous across  = ℎ ! . (BC2-3) In the low layer, the solution of Eq. ( 11) is From (BC1), the constant C 0 is the isotopologue concentration ratio in equilibrium with the liquid sea surface at the sea surface temperature (SST) (Horita et al., 2008), which we obtain from the specified isotopic composition of ocean water and the fractionation factors between liquid water and vapor (Majoube, 1971).Kinetic fractionation is caused by vertically distributed molecular processes concentrated mostly between the surface and z=z * , and is explicitly included by the presence of  ! in Eq. ( 16).This treatment of kinetic fractionation, alone, distinguishes between this IMBL and most other models of 20 atmospheric vapor isotopes near the sea surface.
The second constant of integration in Eq. ( 16) is  !, which is the value of   at z=h 1 .This constant cannot be evaluated at this point, but we return to it shortly.The solutions given by S1 and S2 are long expressions that are far less amenable to evaluation and interpretation than Eq. ( 16), their equivalent for the low layer.Furthermore, they still contain constants  ! and  !, introduced from Eq. ( 16) via the BC's for continuity at z=h 1 .Thus, the solutions for the middle and high layers cannot be evaluated until after  ! has been 5 found.
In order to find  !, it is necessary to apply (BC6) to equations in S1 and S2 (Supporting Information).The somewhat lengthy result is the equation in S3 in Supporting Information.Once  ! has been evaluated, it is feasible (but tedious and slow) to evaluate equations given in S1 and S2 along with Eq. ( 16), completing the evaluation of the unique solution set.

Hybrid Analytical/Numerical Solutions 10
It is more convenient to use a hybrid analytical/numerical approach to finding the solution set.The simple analytic solution for the low layer (Eq.16) can be evaluated in conjunction with a numerical solution for the middle and high layer equations (Eqs.14 and 15).
Numerical boundary value problem solvers normally require the specification of boundary conditions containing only the variables, their derivatives, and numerical constants.Such a solver would not be of use, here, because the constant  ! is not 15 known a priori, so (BC2) and (BC3) cannot be invoked.However, Matlab's © boundary value problem solver bvp5c offers the option of specifying one unknown "parameter" together with two second order boundary value problems and five (instead of the usual four) boundary conditions, and solving for the unknown parameter as well as the continuous variables.
In the analytic problem, this would be equivalent to using 5 boundary conditions to solve for four unknown constants of integration and one unknown "parameter" (C 1 ), essentially what was described in Section 3.1.20 The Matlab © function PBL_analy_numer, in Text S4 of Supporting Information, uses this technique to solve for the isotopologue profiles in the MBL.It calls the solver bvp5c (line 143).The solver bvp5c, in turn, calls the function res (line 416), for the boundary conditions.Since  !appears in (BC2), it can be designated by res as an "unknown parameter", and the five other boundary conditions (BC2-6) can be specified.The boundary value problem that governs the isotopologue profiles in the MBL is thus completely determined.25

Summary
Table 2 contains a list of the eight model parameters that must be specified based on environmental information, and the eight model outputs (either prognostic or diagnostic variables) that are routinely calculated by the model (others can be added).Remember that C is a vector of dimension 3, corresponding to three isotopologues.

Data for Model Validation 5
We use seven published data sets for verification of and comparison with our model output.All of these data sets were collected by shipboard measurements.The summary information is included in Table 3, and cruise tracks are illustrated in Figure 2. Samples from these cruises cover a wide range of the world oceans, from the Arctic Ocean to the northern coast of Antarctica.For earlier data sets, i.e., those by Uemura et al. (2008) and Kurita (2011), samples were collected by the cold trap method, and each sample represents an average of 2-12 hours of vapor trapping.Data from the latter five cruises 10 reported by Benetti et al. (2017) were collected by isotope vapor analyzers with the reported instrument model included in Table 3. Benetti et al. (2017) published data sets with either 15 min or 6 hr resolutions; the 6-hour average data are used for this work.The sea surface temperature (SST), which was either directly measured or estimated by the authors, is reported in all datasets.The relative humidity with respect to SST, RH SST , is either reported (Benetti, et al., 2017) or can be calculated based on the measured air temperature and relative humidity at the sampling height.Both SST and RH SST are important 15 variables in our model validations.

Distribution of Parameters for Verification Runs
In this section, we discuss the ranges of parameter values used in the IMBL model verification simulations.The values are summarized in Table 4.

Sea Surface Temperature (SST) 20
The range of SST used in the simulations was from -2 to +30℃, covering the range of the cruise data sets in Table 3.

Heights h 1 , 2 , 3
A finite span of values was not used for either h 2 or h 3 , because results are insensitive to both, and computations were thus reduced in number.The single value used for the MBL height (h 3 ) was 1000 m, a typical MBL height (Stull, 1988), especially in convergent vapor source areas.Similarly, 650 m was the only value used for h 2 .On the other hand, a full range 25 of values was used for h 1 , because an informal survey of marine radiosonde data suggests that h 1 may range from 50 to 200 m and our results are sensitive to the value of h 1 .

Eddy Diffusivity K max
The eddy diffusivity, K, in the atmosphere boundary layer varies widely over many orders of magnitude.Stull (1988) cited values from 0.1 to 2000 m 2 s -1 , with typical values on the order of 1 to 10 m 2 s -1 for the atmosphere boundary layer.Olivié et al. (2004) presented a calculated range of 0.01 to 3000 m 2 s -1 in the lowest 3 km of the atmosphere for 15 days in July, 1993 at two continental and one marine locations; their maximum value (K max ) above the Pacific ocean location ranged from about 5 3 to 300 m 2 s -1 .Holtslag and Boville (1993) reported calculated zonal and 31-day average eddy diffusivities between 60°S and 60°N; K max ranged from 20 to 60 m 2 s -1 .For K max greater than 10 m 2 s -1 , model isotopic ratios change only negligibly.At K max values less than 1 m 2 s -1 , the kinetic isotopic fractionation increases significantly as K max decreases.We, therefore, use a K max range from 0.01 to 100 m 2 s -1 , to obtain the full extent of kinetic fractionation.

Properties of Subsiding Air (r E , C E , β) 10
Modeling of convergence requires knowledge of the mixing ratio of the descending air (r E in g vapor kg -1 dry air) and its isotopic compositions (C E ), as well as its proportion (β) in the air converged into the MBL.Recall that C E is a vector of length three, corresponding to the concentration ratios of the three modeled isotopologues.The C E value of H 2 16 O is only very slightly less than r E , while values of C E for the other two isotopologues can be obtained from r E and isotopic ratios (δD and δ 18 O) of the vapor.15 Vertical profiles of r E over the ocean have been well observed.We used standard resolution radiosonde data from the University of Wyoming (http://weather.uwyo.edu/upperair/sounding.html) to examine typical tropospheric values and vertical profiles of the mixing ratio.Generally, the mixing ratio decreases rapidly with height within the lower troposphere, and becomes quite small above the mid troposphere.Subsiding air originating in the high troposphere has a correspondingly low mixing ratio.For example, at 500 hPa, the summer-averaged mixing ratio value in the ECMWF (European Centre for 20 Medium-Range Weather Forecasts) data varies from 0.5 to 2 g kg -1 .Most cruise data in Table 3 were obtained between summer and fall, particularly the high latitude ones, and we thus use a range of 0.5 -2 g kg -1 for r E (Table 4) in the simulations.
Measurements of the isotopic composition of vapor are scarce at high altitude.Worden et al. (2007) determined the isotopic composition of tropospheric water vapor from global satellite observations.Values of δD averaged over the altitude range 25 corresponding to pressures between 800 and 550 hPa were found to vary from -180‰ to -250‰ over the extra-tropical ocean.A more recent update reported δD values from -140‰ to -250‰ between 900 and 425 hPa (TES v5 from Sutanto et al., 2015).Ehhalt's (1974) measurements from aircraft above the Pacific Ocean offshore of Santa Barbara, California showed vertical variations of δD from -96‰ to -462‰ between 15 to ~10,000 m for all seasons.The averages for all seasons range from -205 at 800 hPa to -303‰ at 550 hPa.Ehhalt's range is lower, but overlaps the range of satellite values 30 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.(Worden et al., 2007;Sutanto et al., 2015).There are no corresponding measurements of δ 18 O.For the verification simulations, we use a representative value of -239‰ for δD, and -33‰ for δ 18 O (Table 4).Although this choice is somewhat arbitrary, we show that it is adequate for most cruise data sets.To demonstrate the effect of this value, we also show model results with δ 18 O of -28‰, as a comparison.
The proportion, β, of mid-tropospheric air within air converged into the modeled column of the MBL varies with 5 atmospheric conditions including MBL stability, wind speed, and surface roughness.We use a range of values for β from 1% to 10%, which are conjectured values, in the verification simulations.

Upward Velocity (w/ρ)
NCEP/NCAR reanalysis data (Kalney et al., 1996) show that the upward velocity at 850 hPa ranges globally from 0.01 to 0.4 Pa s -1 in magnitude.A typical value of w is 0.1 Pa s -1 over the ocean in summer.We use a range from 0.012 to 0.18 Pa s -1 , 10 corresponding to 0.01 to 0.15 m s -1 .

Other Parameters and Constants
In addition to the parameters discussed above, a few more parameters and/or constants are needed for the simulations.For the isotopic compositions of ocean water, both δD and δ 18 O are set to zero.The molecular diffusivity of H 2 16 O in air is assumed to equal that of bulk water vapor, whose temperature dependence in m 2 s -1 is given by the polynomial fit to Bolz and 15 Tuve's (1976) data (Nellis and Klein, 2009) Merlivat (1978).The turbulent diffusivity at the top of the MBL is set to 100 K m ; while there is little data with which to justify this choice, it suffices because the results are insensitive to it.
The values listed in Table 4 yield 2835 combinations, the result of which is the set of model results we discuss in the next 20 section.

Results and Discussion
In this section, we discuss the characteristics of the model output and their physical significance, and compare the output with observations.We first show vertical profiles of isotopic properties of vapor in the MBL for a representative set of parameters, and then we present the entire set of result of 2835 calculations.These results are then compared with cruise data 25 in both δD vs. δ 18 O and d-excess vs. δ 18 O spaces.We end the section with a discussion of model limitations and potential future developments.

Characteristics of Model Results and Model Validation
While a full discussion of parameter sensitivities and the associated physical processes is the subject of an anticipated companion paper, we point out a few major features of the model output that will guide our discussion of model validation.

Vertical profiles
As a one-dimensional model, the IMBL model yields the vertical distribution of the isotopic quantities δD, δ 18 O, and d-5 excess (=δD-8δ 18 O).Figure 3 illustrates a typical result.Vapor isotopic values δ 18 O and δD are both high near the sea surface, where vapor is in equilibrium with ocean water.With increasing height, isotopic ratios and humidity decrease because of the mixing of MBL vapor with isotopically depleted vapor that descends from the upper atmosphere outside, and then is converged into, the modeled column.The upper atmosphere vapor has much lower values of both δD (-239‰) and δ 18 O (-33‰), but a higher value of d-excess (25‰), than vapor in equilibrium with ocean water.10 The profiles in Figure 3 display strong curvature with very steep gradients near the sea surface, diminishing to negligibly small gradients throughout the MBL.This curvature arises from the rapid change of K from very small molecular values within the thin laminar layer near the water-air interface to large turbulent values above the laminar layer.In this work, the thickness of this layer is characterized by z * , the height of the crossover between molecular and turbulent diffusivities, below which turbulent diffusion is suppressed (See Eqs. 4 and 5).15 The small molecular diffusivity that dominates diffusion in the laminar layer --in particular, its differences among isotopologues --is the cause of kinetic fractionation.Kinetic isotope fractionation is reflected by d-excess that changes more sharply near the surface than does either δD or δ 18 O.The smaller inset of d-excess vs. height plot shows variations within 20 cm of the water-air interface.The z * value, which is 2.7 cm in this particular run, is indicated in the inset by the dashed line.
The effect of turbulent diffusion increases with height, and thus the rate of change in d-excess with height decreases rapidly 20 as the height increases.
Most in situ observations are conducted at a fixed height above the sea surface.The seven cruise data sets (Table 3) were collected at heights between 10 and 20 m.In these cases, each measurement represents an air sample at a given height along a vertical profile.As shown in the one calculation depicted in Figure 3, isotopic gradients are greatest near the sea surface; in this example, over just 15 m (which is only 1.5% of the total height of the MBL) δ  3).In the upper small inset, superimposed in red over the 15 m values are isotopic ratios at both 0 5 m (in equilibrium with seawater at 5 °C) and 15 m for the particular simulation presented in Figure 3, giving a different perspective on the vertical isotopic change.Vapor at 15 m for this particular run has about average deviation from the sea surface equilibrium vapor.Other runs may have larger or smaller vertical gradients in either or both δD or δ 18 O.As discussed earlier, most changes occur below 10 m; above 10 m the change in isotopic composition is relatively minor.
The lower small inset in Figure 4 shows a comparison of two sets of simulations (2835 runs each) using different oxygen 10 isotopic ratios for the upper atmosphere air.Only the boundaries of the output areas are shown, with blue being identical to the main graph, and red indicating the range of results produced using -28‰ (rather than -33‰) for the δ 18 O value of the upper atmosphere vapor.
The output in δD-δ 18 O space (Figure 4) defines a quadrilateral with each corner labeled A through D. The edges (BC, CD, DA and AB) have specific physical significance.Line BC (line b) connects all points (squares) representing isotopic values 15 of vapor in equilibrium with seawater, for the range of sea surface temperatures considered.With increasing sea surface temperature, the points shift from lower left (C) to upper right (B).Points close to this line reflect model parameters that permit very little kinetic isotopic fractionation to occur between the sea surface and 15 m, and very little influence of descending air (whose isotopic composition is point E).Close examination reveals that the points near line BC were generated with the largest turbulent mixing coefficients (highest K max ), and a very small fraction of external air (small 20 β~0.01).Consequently, z * values are very small (~1x10 -5 m), and the relative humidity with respect to SST, RH SST , is close to saturation, both of which are responsible for the small kinetic effect.Line CD bounds points that have the smallest deviation from the straight line CE (line c) that represents mixing of vapor in equilibrium with SST at the coldest temperature considered (-2 o C , point C) and vapor from the descending high altitude dry air (E).Increasing contribution from kinetic isotopic fractionation moves points increasingly above this line (see further 25 discussion below).Therefore, points on this line represent no kinetic fractionation, with the influence of upper atmosphere air increasing from C to E. In other words, if the SST is -2 o C, line CE represents a lower bound on isotopic mixing.At a fixed SST and ocean isotopic ratio, this line rotates with changing isotopic ratios in the air aloft, for example, line CF in the lower inset of example, if the SST is 30 o C, then the mixing line would be BE (not shown), and all isotopic ratios of vapor evaporated from this sea surface should be above this line.
The points along line AB represent vapor evaporated from ocean water at SST=30 o C. Their spread reflects the influence of kinetic fractionation; moreover, they are not significantly influenced by mixing with upper atmosphere air.This inference is supported by their small values of K max (0.1 m 2 s -1 ), large z * (0.3-0.5 m), and low r E of the upper atmosphere (0.5 g kg -1 ). 5 Theoretically, the slope of AB should have a limit of 0.88 (shown as line a), the ratio of the kinetic fractionation factors of HDO and H 2 18 O (25.1 and 28.5‰, respectively, because K m /K * m =1.0251 and 1.0285, respectively, where the star represents the heavy isotopologue; Merlivat, 1978).With the set of parameters in Table 4, the slope of AB is about 1.5, slightly greater than its lower limit (0.88).Therefore, line a sets the upper bound vapor isotopic ratios for SST of 30 o C. In other words, the theoretical limit for the highest isotopic ratios at a given SST should be along a line that starts from a point representing 10 vapor in equilibrium with seawater (δ 18 O=0, δD=0 in this calculation) at this temperature and extends to the lower left with a slope no less than 0.88.
Line AD bounds isotopic compositions reflecting the entire range of SST values; both kinetic fractionation and mixing with the upper atmosphere have significant influences on these points.The ambient conditions are characterized by small K max (0.01 m 2 s -1 ), large β (0.1) and relatively high z * values (0.1-0.5 m).This AD "line" is not as strictly defined as other lines in 15 that it does not have an absolute theoretical limit and so may change with the range of parameter space.In subsequent discussion, we refer to line a as the upper limit, line b the side limit, and line c the bottom limit of the vapor distribution in the δD-δ 18 O space, consistent with their positions in Figure 4.
In summary, the shape of the output in Figure 4 is controlled by three factors, 1) the SST, 2) the degree of kinetic isotopic fractionation, and 3) the influence of upper atmosphere air.While SST is relatively independent of other factors, kinetic 20 fractionation and effect of upper air depend on various combinations of atmospheric conditions, including the intensity of turbulent mixing (K max ), the mixing ratio of the descending air (r E ) and its isotopic ratios, the proportion of upper atmospheric air advected into the evaporation column (β), and the vertical velocity (w).Note that in this model, the relative humidity with respect to SST (RH SST ) is not, and cannot be, prescribed.On the contrary, it is an outcome of the same meteorological conditions of the MBL that affect the isotopic ratios, although it also feeds back on kinetic isotopic 25 fractionation by controlling the vertical gradient for vapor diffusion.
Model output and observational data for each individual cruise are compared in Figure 5. Model output is calculated at the observation height of the corresponding cruise, indicated in the graph.Also included in each plot are compositions of vapor in equilibrium with seawater at the lowest and highest SSTs of the cruise.The theoretical borders under specific cruise conditions are shown as solid lines; observed isotopic ratios are expected to fall within these theoretical limits (if consistent 30 with the assumed ocean water and descending air isotopic ratios).We make the following observations of Figure 5. First, the vast majority of the observed data (~95%) do indeed fall within the expected range.This confirms not only the successful conceptualization of the model but also that our choices of parameter values are reasonable.
Second, in all seven data sets, the influence of the isotopically depleted vapor from descending air is demonstrated by points with low isotopic ratios.These compositions are difficult to model using C-G type models, particularly with the invalid 5closure assumption (e.g., Jouzel and Koster, 1996;Benetti et al. 2015).This result highlights the importance of convergence in affecting boundary layer vapor isotopic ratios, as it introduces dry, depleted air from aloft into the MBL.Such influence of upper atmosphere air on the boundary layer has been recognized by Benetti et al., (2015Benetti et al., ( , 2018)), although for quiescent subsidence regions that our model does not treat.Furthermore, their conceptual treatment is still based on the closure assumption.10Third, for the ACTIV cruise (Figure 5c), a number of points fall below the lower limit, suggesting that the isotopic ratios of the descending air we used for the simulation may not be representative in area during the observation period.The mismatch suggests a value that is more enriched in 18 O, or depleted in deuterium, or both, than the values used for the simulation.
Fourth, in four cruises (PIRATA, STRASSE, Bermuda and RARA; Figure 5d-g), there are points that are above the upper 15 limit.However, in all cases except RARA, the enrichment above the upper limit is small in magnitude, and may be explained by slight variations in seawater isotopic ratios.For RARA, however, the enrichment above the upper limit is significant.One possible explanation is the influence of sea spray.When describing the sampling conditions, Benetti et al. (2017) particularly noted that they could not completely rule out the contribution of small droplets of sea spray to the vapor composition.
However, such an influence seems relatively small considering the great leverage of seawater isotopic composition.Figure 20 5h shows the direction and magnitude of sea spray influence; the mixing of sea spray droplets should cause enrichment such that the data would be distributed in the triangular area bordered by the dashed lines.Detailed examination of cruise logs in the future will be helpful to confirm and quantify the sea spray contribution to MBL vapor.In summary, by comparing calculated values and observational data in δD-δ 18 O space, we conclude that the model is remarkably successful.We pointed out three factors that may cause observations to fall outside the predicted range, namely 25 1) variation in ocean water isotopic ratios, 2) variation in the isotopic ratios of the upper atmospheric vapor, and 3) influence of sea spray on vapor isotopes.In section 6.2, we discuss several other model assumptions that may limit the consistency between model results and observations.Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.

Deuterium excess (d-excess)
The relationships between d-excess and both sea surface temperature (SST) and relative humidity with respect to SST (RH SST ) have been intensively discussed by the isotope hydrology community.Originally defined by Dansgaard (1964) for precipitation as δD-8δ 18 O, d-excess has been used to infer conditions at the moisture source location.It has been demonstrated that d-excess varies with SST and inversely with RH SST (Johnsen et al, 1989;Petit et al., 1991).Later 5 investigators have used these concepts to infer changes in moisture source conditions recorded in polar ice (e.g., Vimeux et al., 1999;Masson-Delmotte et al., 2005a, 2005b).
The relationships between d-excess and SST, and between d-excess and RH SST , are shown in Figure 6.Our model, as expected, exhibits a significant dependence of d-excess on both SST and RH SST .Regression of d-excess against SST explains 16% of the variance in d-excess, with a coefficient of 0.35‰ o C -1 .Regression against RH SST explains 78% of the 10 variance in d-excess, with a coefficient of -0.43‰ % -1 .These values are very similar to d-excess sensitivities of 0.35‰ o C -1 to SST and -0.45‰ % -1 to RH SST , respectively, cited by Vimeux et al., (1999) based on calculations by Johnsen et al., (1989).
All three processes discussed earlier, i.e., changing SST, degree of kinetic fractionation, and extent of mixing with the subsiding air, result in changes in d-excess.This is seen by the fact that the theoretical lines in Figure 6 representing each of 15 the three processes have non-zero slopes.Although at the sea surface d-excess increases with SST, a much larger spread occurs at 15 m due to the height-dependent influence of the descending air and kinetic fractionation.For each value of SST, the points at 15 m form a triangular area, within which d-excess varies significantly (Figure 6a).Such a triangular distribution of isotopic data in the d-excess vs. δ 18 O space has been reported by Steen-Larsen et al. (2015) for observations off the coast of Iceland.This two dimensional distribution explains the significant, though relatively poor, correlation 20 between d-excess and SST.
Figure 6b shows that the lowest RH SST tends to correspond to the highest d-excess distributed near corner A. However, near point D, where d-excess is also relatively high compared with values at the water-air interface, RH SST is relatively high and kinetic fractionation is limited.Therefore, while d-excess tends to increase as RH SST decreases, the relationship is not one to one (note how color changes horizontally in Figure 6b).Another way to see this is to trace the color change along lines 25 parallel to CB in Figure 6a for changes in SST, and parallel to AB in Figure 6b  outside the predicted region were discussed earlier, and we do not repeat that discussion here.Second, the dependence of the observed d-excess on RH SST and SST, as shown by the color distributions, is very similar to that of model calculations.For SST, the color along lines parallel to CB changes from green to red with a d-excess increase.Similarly, RH SST values are relatively high near lines CB and CE, and decrease (with significant noise) towards corner A as d-excess increases.Finally, as predicted, the observed d-excess correlates significantly (p<0.0001) with SST and with RH SST .The sensitivity of d-excess 10 observations to RH SST is -0.36‰ % -1 , comparing favorably with corresponding model sensitivity of -0.43‰ % -1 .For SST, the sensitivity for observations is 0.15‰ o C -1 , about half of that predicted by the IMBL model (0.35 ‰ o C -1 ) using SSTs ranging from -2 to 30 o C. Detailed discussion of the sensitivity differences between simulations and observations is beyond the scope of this paper, and a full understanding of these relationships should be an important goal for future investigations.

Model Applicability, Limitations and Future Development 15
The IMBL model described here has considerable potential value for many isotope hydrology applications.First, as vapor isotopic measurements become increasingly available, application of the model at different locations and times of year provides a vehicle to explore and understand relationships between meteorological conditions and isotopic compositions of vapor both within the MBL and delivered to the free atmosphere.Comparisons of simulations and observations with much more intensive observations than cited in this work may be conducted.For example, during an isotopic vapor measurement 20 campaign, measured variations of the isotopic composition of ocean water and vapor may be combined with model calculations to constrain the diffusion coefficient, which may then be related to sea surface roughness, wind speed, vertical velocity, and sea spray occurrence.Second, the output of this model, i.e., the isotopic ratio of vapor delivered to the free atmosphere, can be used to provide initial conditions for Rayleigh-type condensation or transport models.The sensitivity of precipitation isotopic ratios to the range of meteorological conditions at the moisture source region and their change over 25 time and space can be investigated for modern hydrological cycles in association with atmospheric circulation, as well as under ancient climate conditions.Third, an understanding of physical processes important for controlling the isotopic compositions of water (both vapor and precipitation), gained from these applications, can be used to improve the physical representation of marine boundary layer processes in GCMs.This IMBL model may not adequately describe several meteorological scenarios, and its use under those conditions should 30 be made with caution.First, the model requires that the air column in the model domain have a positive (upward) vertical velocity, i.e., air must be converging and rising.This assumption is made to ensure that the model column represents a Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.moisture source area, delivering vapor to the free troposphere where it will ultimately produce precipitation.If the vertical velocity is negative (subsiding air), the air in the column diverges rather than converges, meaning that the mass conservation equations used here would not be correct.Either way, the important outcome is that upper atmosphere vapor is incorporated into the MBL.It is possible that isotopic distribution changes somewhat with specific mixing scenario, while the theoretical limits remain the same as shown by this work.Second, the model does not include exchange between vapor and liquid if air 5 is supersaturated, forming clouds or precipitating.Third, as discussed earlier, the model does not include the effects of sea spray.Fourth, the modeled column is not subject to horizontal advection (except for convergence).Fifth, the IMBL is a steady state model.
We envision future improvements in the IMBL model.In particular, we anticipate generalization of the model to include areas of divergence (descending air), areas with sea spray, and/or terrestrial areas.10

Conclusions
We have constructed a new isotope marine boundary layer (IMBL) model to calculate the isotopic composition of vapor in the marine boundary layer as well as that of vapor lifted to the free atmosphere.The model divides the MBL into three layers, each with its own transport characteristics.Compared with earlier Craig and Gordon (1965) type bulk evaporation models, this 1-D (vertical) model makes the following improvements: 1) It explicitly includes the process of horizontal 15 convergence in the middle layer; convergence provides mass to balance the upward advection supplying moisture for cloud formation and precipitation in the free atmosphere.This formulation requires specification of the properties of subsiding air that mixes with low altitude air and converges into the model column.2) The eddy diffusion coefficient is height-dependent, equal to the molecular diffusion coefficients for each isotopologues at the air-water interface, and increasing linearly through the lower layer to a maximum value, K max , remaining constant through the middle layer, and decreasing linearly to K t over 20 the top layer.The advantages gained from these improvements include 1) the model solves for both isotopologue concentrations in and fluxes through the MBL, not just fluxes; 2) kinetic isotopic fractionation becomes a diagnostic variable rather than a required parameter, without adding to the total number of parameters (degrees of freedom) of the model; 3) relative humidity is also no longer a specified parameter, but rather becomes a diagnostic variable, thus enabling the use of the model to investigate and possibly predict evaporation rates; 4) calculation of vertical profiles of concentrations and 25 fluxes of isotopologues (or isotopic ratios), providing an opportunity to compare model output with observations at a specific height or multiple heights; and 5) the air at the top of the MBL (at z = h 3 ) is the air mass supplied to the beginning of a Rayleigh trajectory, which can be used in many isotope studies.troposphere.Kinetic fractionation near the sea surface causes d-excess to increase rapidly with height, particularly between the air-sea interface and height z * , where molecular diffusion dominates over turbulent mixing.
Model simulations using reasonable ranges of parameters are validated using seven sets of shipboard MBL observations.The resulting range of δD and δ 18 O forms a quadrilateral-shaped pattern in the δD-δ 18 O plane.Three processes generate boundaries for the quadrilateral, or constraints on the isotopic ratio distributions, including 1) the set of vapor isotopic ratios 5 in equilibrium with seawater at various SST 's (right side boundary), 2) mixing between vapor descended from the upper atmosphere and vapor in equilibrium with seawater at the air-water interface (lower boundary), and 3) kinetic isotopic fractionation that starts with equilibrium vapor and extends to more depleted values of δD and δ 18 O, with a slope no less than 0.88 (upper boundary).About 95% of observations fall into the theoretically predicted quadrilateral region, demonstrating the success of the model 10 conceptualization and parameter choices.This remarkable agreement highlights the importance of convergence and entrainment of descending, isotopically depleted air to boundary layer isotopic ratios.This feature is new to this model, and was not considered in earlier simple models, although some (e.g., Benetti, 2015) do include mixing by mathematically unspecified mechanisms other than convergence, in meteorologically regions distinct from those we consider.The simulation-observation comparisons also point to at least three factors that may explain the 5% of data falling outside the 15 predicted region, including 1) variations in seawater isotopic ratios, 2) inaccurate choice of isotopic ratios for the subsiding air, and 3) influence of sea spray.It is also possible that meteorological scenarios not explicitly considered by the model are responsible for the relatively minor model-data mismatch.Such factors may include low level air divergence (downward vertical velocity in the middle and upper MBL), vapor-liquid exchange (during precipitation events or within clouds), and the presence of lateral advection.20 Simulated d-excess compares remarkably well with observations.While the effects of sea surface temperature (SST) and relative humidity with respect to SST (RH SST ) are well-understood, we draw attention to the influence of the upper atmosphere air in controlling d-excess in the marine boundary layer.In this simulation, the d-excess value of the descending air is greater than that of vapor in equilibrium with seawater, and the contribution of the former to MBL air results in an increase in d-excess of vapor, even in the absence of kinetic isotopic fractionation.25 The IMBL model can be used in a number of ways.First, numerical experiments with the model help to better understand the effects of boundary layer processes and climatic conditions on isotopic compositions of vapor within and vapor fluxes through the MBL.A second application could be to investigate how temporal and spatial differences in moisture source regions affect the isotopic composition of remote precipitation under both modern as well as paleo-climate conditions.Third, it is important to investigate the relationship between MBL isotopes and evaporation rate and, perhaps, to develop methods 30 to measure the latter indirectly from simultaneous observations of isotopes and meteorological conditions.Finally, the Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
The model introduced here is a one-dimensional (vertical) model with three layers within the MBL.It adopts the following enhancements to improve upon the earlier, classical models: 1) It explicitly includes vertical velocity and horizontal convergence of air and vapor, notwithstanding the difficulties of specifying the fluxes and isotopic properties of converging air.2) It uses a height-dependent eddy diffusion coefficient without increasing the total number of free parameters (degrees 30 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License. of freedom) in the model.3) It does not make the closure assumption that isotopic flux equals isotopic composition.4) It solves not only for isotopologue fluxes, but also concentrations.5) MBL humidity and kinetic fractionation factors are no longer required input parameters but are calculated.6) Vapor fluxes are no longer constant with height.Above, we have made several references to applying an MBL model to the initiation of a Rayleigh model of vapor trajectories in the free troposphere, but there is another crucial role for an MBL model.It is the model's application to 5 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.Sketched on the right side of Figure 1 are vertical profiles of the diffusion coefficient K(z) and the (dynamic) vertical velocity w(z).The profile of K(z) is consistent with typical variation with height of the eddy viscosity diffusion coefficient in the boundary layer, based on O' Conservation of mass for an isotopologue affected only by diffusion, temporarily neglecting convergence and advection, takes the form Atmos. Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
Similar to the low layer, the middle and high layers have analytic solutions.As is standard with boundary condition problems, the general solutions are found first.Then BC's 2-5 are introduced into the solutions, and the four new constants Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License. of integration are solved for (in terms of the model parameters).The results are shown as Texts S1 and S2 in Supporting Information, respectively.
Figure 4 is somewhat arbitrary, but is approximately the average of the observation heights, that range from 10 to 20 m, for the seven data sets (Table3).In the upper small inset, superimposed in red over the 15 m values are isotopic ratios at both 0 5 Figure 4. Similarly, mixing lines between equilibrium vapor at higher SST's should be straight lines connecting point E and points along line b representing vapor in equilibrium with seawater at different temperatures.For 30 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
for changes in RH SST .Interestingly, along CD, neither SST nor RH SST varies significantly, regardless of substantial variation in d-excess.Obviously, RH SST is not the sole influence on d-excess, and even the combination of both RH SST and SST does not completely determine d-excess in the MBL.Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
The standard output of the model includes vertical profiles of δD, δ 18 O and d-excess.Near the sea surface, δD and δ 18 O are high and close to equilibrium with the ocean water, and d-excess is low.With increasing altitude, δD and δ 18 O decrease due 30 to both kinetic fractionation and mixing with converging isotopically depleted air that contains subsided air from the free Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-709Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 11 September c Author(s) 2018.CC BY 4.0 License.

Figure 1 :
Figure 1: Cartoon of the MBL model.The modeled region is indicated by the cylinder.It is comprised of three layers -the low von Kármán layer, the middle convergence layer, and the top stabilizing layer.The heavy straight arrows represent the flow of semi-moist air ascending through the middle and top layers, and through the top of the MBL into the free atmosphere (above the 5

Figure 2 :
Figure 2: Map showing the tracks for seven cruises that generated vapor isotopic observations used in this work.Information about each cruise is listed in Table3and the associated references.

Figure 3 :
Figure 3: Model simulation of vertical distributions of δ 18 O (left), δD (middle) and d-excess (right).Parameters are SST=5 o C, K max =0.1 m 2 s -1 , h 1 =120 m, r E =0.5 g kg -1 , w a =0.15 m s -1 , and b=0.05.The horizontal solid lines mark the heights of h 1 and h 2 (120 and 650 m).The inset graph shows d-excess variation with height in the 20 cm just above the sea surface.The dashed line marks 5 the value of z * (0.027 m), below which molecular diffusion is more significant than turbulent diffusion.

Figure 4 :
Figure 4: δD vs. δ 18 O relationship at 15 m height for 2835 model calculations (+).The output defines a quadrilateral with corners labeled by A-D.Also shown are δD and δ 18 O values of the descending air (E, u), and isotope values of vapor in equilibrium with the seawater (o ,along line b) at SSTs of -2, 5, 10, 15, 20, 25, and 30 o C from C to B, respectively.Solid lines labeled a, b, and c bound the theoretical limits of vapor isotopic values, with b being vapor in equilibrium with seawater, c being a mixture of vapor 5

Figure 5 : 5 highest
Figure 5: Comparisons between simulated and observed isotopic ratios for each of the seven cruises (a-g), and a redraw of the RARA data with ocean water composition indicated (h).For each cruise, the simulated values are calculated at the height of the observations, indicated in the plot.If not otherwise indicated, calculated isotopic values are shown as crosses (+), the values of the descending air as a diamond (u), and observations as red circles (¢).For each cruise, vapor in equilibrium with lowest and 5

Figure 6 : 5 6
Figure 6: Plots of d-excess vs. δ 18 O, showing relationship between d-excess and SST for the simulation (a) and model-data comparison (c), and between d-excess and RH SST for simulation (b) and model-data comparison (d).As in Figure 4, isotopic values are calculated at 15 m height.Points corresponding to those from A to E in Figure 4 are also shown, with point E being the isotopic composition of the descending air.Straight lines are theoretical limits for processes labeled in (a) and (b) (also discussed in section 5