Isotopic modeling of the sub-cloud evaporation effect in precipitation

Science of the Total Environment 544 (2016) 1059–1072

Contents lists available at ScienceDirect

Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv

Isotopic modeling of the sub-cloud evaporation effect in precipitation V. Salamalikis a,⁎, A.A. Argiriou a, E. Dotsika b a b

Laboratory of Atmospheric Physics, Department of Physics, University of Patras, GR 26500 Patras, Greece Stable Isotope Unit, Institute of Nanoscience and Nanotechnology, National Center of Scientific Research ‘Demokritos’, Ag. Paraskevi Attikis, 15310 Athens, Greece

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T

• Modeling the sub-cloud evaporation effect in terms of the isotopic composition. • Two isotopic modeling approaches used a) isotope-mixing evaporation model and b) numerical isotope evaporation model. • Model evaluation using data from Madrid and Vienna GNIP stations. • Vertical distribution of stable isotopes in falling raindrops in terms of isotopically transformed difussive equations.

a r t i c l e

i n f o

Article history: Received 17 September 2015 Received in revised form 9 November 2015 Accepted 14 November 2015 Available online xxxx Editor: D. Barcelo Keywords: Stable isotopes Falling raindrops sub-cloud evaporation effect Isotope evaporation models

a b s t r a c t In dry and warm environments sub-cloud evaporation influences the falling raindrops modifying their final stable isotopic content. During their descent from the cloud base towards the ground surface, through the unsaturated atmosphere, hydrometeors are subjected to evaporation whereas the kinetic fractionation results to less depleted or enriched isotopic signatures compared to the initial isotopic composition of the raindrops at cloud base. Nowadays the development of Generalized Climate Models (GCMs) that include isotopic content calculation modules are of great interest for the isotopic tracing of the global hydrological cycle. Therefore the accurate description of the underlying processes affecting stable isotopic content can improve the performance of isoGCMs. The aim of this study is to model the sub-cloud evaporation effect using a) mixing and b) numerical isotope evaporation models. The isotope-mixing evaporation model simulates the isotopic enrichment (difference between the ground and the cloud base isotopic composition of raindrops) in terms of raindrop size, ambient temperature and relative humidity (RH) at ground level. The isotopic enrichment (Δδ) varies linearly with the evaporated raindrops mass fraction of the raindrop resulting to higher values at drier atmospheres and for smaller raindrops. The relationship between Δδ and RH is described by a ‘heat capacity’ model providing high correlation coefficients for both isotopes (R2 N 80%) indicating that RH is an ideal indicator of the sub-cloud evaporation effect. Vertical distribution of stable isotopes in falling raindrops is also investigated using a numerical isotopeevaporation model. Temperature and humidity dependence of the vertical isotopic variation is clearly described by the numerical isotopic model showing an increase in the isotopic values with increasing temperature and decreasing RH. At an almost saturated atmosphere (RH = 95%) sub-cloud evaporation is negligible and the isotopic

⁎ Corresponding author. E-mail address: [email protected] (V. Salamalikis).

http://dx.doi.org/10.1016/j.scitotenv.2015.11.072 0048-9697/© 2015 Elsevier B.V. All rights reserved.

1060

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

composition hardly changes even at high temperatures while at drier and warm conditions the enrichment of 18Ο reaches up to 20‰, depending on the raindrop size and the initial meteorological conditions. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Precipitation is the most important component of the hydrological cycle for the determination of the water budget over a certain area. The accurate description of the precipitation patterns as well as the underlying mechanisms that modulate precipitation is a crucial issue for the robust parameterization of climate models and to further improve the quality of the precipitation forecasts. In arid and semi-arid climates sub-cloud evaporation affects the falling hydrometeors during their descent from the cloud base towards the ground modifying the size of the raindrops reaching the ground. The use of environmental tracers sensitive to such changes such as stable water isotopes (18O and 2H) can contribute to the better understanding of the presence and the impact of the sub-cloud evaporation effect in precipitation. For example, Worden et al. (2007), using satellite data, reported that 20% (and up to 50%) of the rainfall is re-evaporated under tropical convective clouds. Stable water isotopes have been widely used as tracers since their values reflect the water transport and the various phase changes occurring at different spatio-temporal scales in the hydrological cycle. Driven by their mass and vapor pressure differences, the physical properties of the different isotopic species vary according to the environmental conditions leading to changes in their abundance due to equilibrium and kinetic fractionation effects (Dansgaard, 1964). Due to the lower vapor pressures, the heavier isotopic species are preferentially accumulated in the condensed phases information on the mechanisms related to cloud and precipitation formations (Jouzel, 1986). The isotopic composition of water is expressed via the delta (δ) notation in per mille (‰) which is the fractional deviation from unity of the isotopic ratio of a water sample (RiS) relative to that of a standard (RiVSMOW); here the Vienna Mean Ocean Water (VSMOW) (Eq. (1)). i

δ ð‰Þ ¼

RiS RiVSMOW

! −1  103 :

The falling raindrops tend to attain equilibrium with the ambient air enriching isotopically the condensate phase of water in comparison to the initial isotopic composition at the cloud base (Stewart, 1975; Lee and Fung, 2008). Isotopic exchange is more pronounced for smaller raindrops. The adjustment/equilibration time - the time a droplet needs to reach the 1/e of the isotopic composition of the environmental vapor - is positively related to the raindrop size indicating that a lower equilibration time and a higher equilibration degree are associated to lower sized raindrops (Friedman et al., 1962; Lee and Fung, 2008). The reduction of the raindrop size due to evaporation depletes the surrounding air at any falling step on the assumption that the isotopic composition of the ambient atmosphere is lower than the isotopic composition of the water vapor near the raindrop's surface. The distribution of stable isotopes in precipitation is modeled via the Rayleigh distillation process. The temperature and rainfall amount effects are included in a certain extent in the Rayleigh formula in terms of the equilibrium fractionation factor and the rainout fraction. The Rayleigh distillation assumes isotopic equilibrium between liquid and vapor at any step of the precipitation process. Also, precipitation is immediately removed from the cloud base without any interaction with the surrounding atmosphere neglecting the sub-cloud kinetic fractionation phenomena (Criss, 1999). Despite the strong limitations, the Rayleigh distillation is widely applied for the determination of the overall isotopic patterns and the isotopic evolution of precipitation/water vapor along pre-described moisture trajectories (Sodemann et al., 2008). The general purpose of this paper is a discussion for the isotopic modification of falling raindrops due to the sub-cloud evaporation effect. In this work this phenomenon is analyzed using three methods, namely: 1. Meteoric Water Lines (MWLs);

ð1Þ

The symbol i stands for 18O and 2H respectively. The falling raindrops are subject to evaporation during their trip from the cloud base towards the surface through the unsaturated atmosphere. The impact of below-cloud evaporation on the isotopic composition of falling raindrops is strongly related to the prevailing meteorological conditions (temperature and humidity) in the sub-cloud layer, the initial raindrop size and the isotopic composition of the ambient moisture. In case of isotopic equilibrium between a raindrop and the atmosphere surrounding it, the isotopic composition of precipitation is derived by the isotopic composition of the surrounding water vapor and the equilibrium fractionation factor for the liquid–vapor transition. The isotopic composition of the ambient moisture decreases with increasing altitude, following the water vapor distribution in the troposphere. At any vertical level, the isotopic composition of precipitation depends on the isotopic composition of the pre-condensed moisture as well as on the isotopic composition of the ambient vapor (Bony et al., 2008). When RH b 100% and the raindrop surface temperature is lower than the ambient air temperature, a heat flux towards the raindrop delivers the latent heat of evaporation and heats the raindrop. Due to kinetic fractionation an isotopic enrichment is detected in the falling raindrops and a reciprocal depletion in the water vapor, in relation to the humidity conditions and the temperature difference between the raindrop and the ambient environment. The evaporation rate of the raindrop is inversely proportional to the raindrop diameter (~1/D2) implying that smaller raindrops are more favorable to evaporation. In case of saturation (RH ≃ 100%), the phenomenon of isotopic exchange occurs.

2. Isotope-evaporation model (Stewart, 1975; Georgakakos and Bras, 1984); and 3. Numerical isotope-evaporation model (Jouzel, 1986; Zhang et al., 1998). The numerical isotope-evaporation model involves the vertical variations of the raindrop diameter and raindrop temperature due to below cloud interactions as well as the vertical isotope profile of the atmospheric moisture. For the subsequent analysis the increase in the raindrop size due to coalescence and collision and the decrease due to breakup, are assumed negligible (Schlesinger et al., 1988). This model is a simplified and computationally ‘light’ version of the complex isotopic cloud-resolving models applied for the representation of the isotopic composition of water vapor and its condensates in convective updrafts (Bolot et al., 2013) and the investigation of the tropical tropopause layer (Blossey et al., 2010). 2. Data Stable isotopic daily data (18O and 2H) for ground water vapor and rainfall were obtained by the GNIP/ISOHIS database for six stations located in the Northern Hemisphere (IAEA/WMO, 2013). Due to gaps in the meteorological data in the GNIP database, missing values of temperature, precipitation amount and relative humidity were filled using the ECA&D database (Klein Tank et al., 2012). Only two stations (MadridRetiro, Spain and Vienna, Austria) have temperature and relative humidity values available in the ECA&D database while the temperature gaps in Rehovot, Israel were filled using the meteorological information from the nearest station in Bet-Dagan, Israel.

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

1061

raindrops have a uniform spherical shape, the isotopic enrichment due to sub-cloud evaporation is (Froehlich et al., 2008),

3. Methods 3.1. Meteoric water lines and isotopic enrichment of precipitation The various fractionation processes in water can be interpreted using the δ2H-δ18O relationship. Globally, the isotopic composition of precipitation (δ18O and δ2Η) lies along the Global Meteoric Water Line (GMWL), δ2H = 8 δ18Ο + 10‰ (Craig, 1961). In smaller spatial scales (areas and stations) and on an event and longer temporal basis, the isotopic composition of the observed rainfall varies due to a variety of factors projecting the water samples along a new straight line, the Local Meteoric Water Line (LMWL), δ2H = a δ18Ο + b, where both the slope and intercept values may differ from those of the GMWL. In case of sub-cloud evaporation, the values of a and b are significantly lower than those of the GMWL implying that the sub-cloud processes can be identified through the examination of the LMWL (Gat, 2005). In order to describe the deviation from equilibrium, the slope of the LMWL is compared against the Rayleigh equilibrium slope αT (Criss, 1999),    aeq;l−v 2 H δ2 H p −103   aT ¼  aeq;l−v 18 OÞ δ18 Op −103

ð2Þ

with aeq,l − v the equilibrium fractionation factor for the liquid–vapor transition (Table A1). The equilibrium fractionation factors for 18O and 2 H are calculated using the Majoube (1971) formulas in terms of the surface air temperature while the δ-values of precipitation are weighted according to the respective precipitation amount. According to Eq. (2), the Rayleigh slope assumes isotopic equilibrium during the formation of rainfall, neglecting the super-saturation conditions inside the cloud and the secondary evaporation phenomena beneath the cloud base. A higher slope difference (Δα = αT − α) indicates stronger evaporative enrichment and deviation from equilibrium. 3.2. Isotope-evaporation model The detection of the sub-cloud evaporation effect is described also by the difference between the isotopic composition of precipitation in the ground level and cloud base respectively. According to Stewart (1975) the isotopic composition in falling raindrops is given by,   i β Ril ¼ γi Riv;cb þ Ril;cb −γ i Riv;cb f r  n 1−aieq;l−v Dva =Diva ð1−RH Þ ð3Þ aieq;l−v RH ; γi ¼ βi ¼  n  n aieq;l−v Dva =Diva ð1−RH Þ aieq;l−v Dva =Diva ð1−RH Þ where Ril is the isotopic ratio of falling raindrops, Rivcband Rilcb are the initial isotopic ratios for the vapor and liquid phases at the cloud base, aieq , l − v is the equilibrium fractionation factor calculated at the cloud base temperature, RH is the fractional relative humidity, Dva/Diva is the ratio between the molecular diffusivities of the abundant and the rare isotope i, n = 0.58 (Gat, 2000) and fr is the remaining raindrop mass fraction after evaporation. The molecular diffusivity ratios for water isotopes are taken to be 1.032 and 1.016 for 18O and 2H respectively (Cappa et al., 2003). The cloud base temperature is calculated through the mean moist adiabatic lapse rate γ¼−

dΤ ¼ 6:5  10−3 K=m dz

ð4Þ

assuming that the cloud base is at about the 850 hPa (~ 1500 m) isobaric level. Starting from Eq. (3) and assuming that the falling

  Δδð‰Þ ¼ δ zg −δi ðzcb Þ ¼ 103  i

mg fr ¼ ¼ mcb



Dg Dcb

1−

3

γi

! 

aieq;l−v

 βi f r −1 ð5Þ

with D the raindrop's diameter. The indexes g and cb refer to the ground surface and the cloud base respectively. For the determination of the final raindrop mass, a station-precipitation model (Georgakakos and Bras, 1984) is applied,    1=3 4Dva esat ðθw Þ esat ðθd Þ zcb − Dg ¼ Dcb 3 − C 1 Rv Tw Tg

ð6Þ

with Rv = 461.51 J.kg− 1.K− 1 the gas constant for water vapor, Dva (m2·sec− 1) the air diffusivity, esat (Pa) the saturation vapor pressure at the wet-bulb temperature, θw and Tw the wet-bulb temperature in °C and K, θd (°C) the dew point temperature, Tg (K) the air temperature at the ground level, zcb (m) the cloud base elevation and C1 = 7 × 105 kg.m− 3.s − 1 for liquid droplets. The wet bulb and the dew point temperatures are calculated in terms of relative humidity and air temperature (Lawrence, 2005; Stull, 2011) while the MagnusTetens formula (Lawrence, 2005) is applied for the computation of the saturation vapor pressure. Even if the Stewart's model describes the isotopic enrichment of falling raindrops, it has important limitations. The model assumes a homogenous sub-cloud layer in terms of temperature and humidity, neglecting their vertical distribution, considering also that the isotopic contribution in the ambient vapor due to the evaporation of the falling raindrops is negligible. However, those limitations are expectable due to the nature of the Stewart's model, since it is a mixing model between the starting and the final isotopic composition. 3.3. Numerical isotope-evaporation model In order to study the temporal evolution of the isotopic composition in falling raindrops, the rate of change of the isotopic ratio is described by the following equation,      i  dRil d mi =m 1 dmi dm 6 dm i dm ¼ ¼ −Ril −R ¼ : l dt dt m dt dt dt πρw D3 dt

ð7Þ

In Eq. (7), the raindrops are assumed to have a uniform spherical 3

shape with diameter D and mass m ¼ πρw6D . The evaporation rate (dm/ dt) and the rate of change of the isotopic mass (dmi/dt) for the falling raindrops are derived using the mass transfer equation by means of the vapor concentration, dm ¼ 2πDð F v Dva Þ½ρva ðT a Þ−ρvr ðT r Þ dt i  h i dm ¼ 2πD F iva Diva ρiva ðT a Þ−ρivr ðT r Þ dt

ð8Þ

with Ta, Tr the ambient temperature and the raindrop surface temperature respectively. The parameters Dva, ρva, ρvr, Fv, Diva, ρiva, ρivr,Fiv are the diffusion coefficient of air, the vapor densities at the ambient atmosphere and at the raindrop surface temperature, and the ventilation coefficients for the normal and the heavy water respectively.

1062

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

Eq. (7) using Eq. (8) can be transformed as, dRil 3ρvr ðT r Þ ¼ ð f v Dva Þ dt ρ"w D2 !   # ρva;sat ðT a Þ ρva;sat ðT a Þ F iv Diva i i i Rva RH −Rvr −Rl RH −1 :  ρvr ðT r Þ ρvr ðT r Þ F v Dva ð9Þ The isotopic composition of the ambient moisture is given by Riva ¼ ρiva ðT a Þ ρva ðT a Þ

while the isotopic composition at the raindrop's surface is Rivr ¼

ρivr ðT r Þ

ρvr ðT r Þ. At the surface of the raindrop, the vapor and liquid phases are con-

sidered in isotopic equilibrium and the isotopic composition of water vapor can be calculated using the isotopic composition of the liquid phase and the equilibrium fractionation factor at the raindrop temperature, Rivr ¼ ai

Ril

eq;l−v

. The quantity

ðT r Þ

F iv Diva F v Dva

Di

dm ¼ 2πDð F v Dva Þ½ρva ðT a Þ−ρvr ðT r Þ dt  2πDð F v Dva Þ eva;sat ðθa Þ esat ðθr Þ : ¼ RH −  Rv Ta Tr

ð12Þ

Using the Eqs. (9) and (11) the vertical change of the raindrop diameter equals  dD 4ð F v Dva Þ eva;sat ðθa Þ esat ðθr Þ : − ¼  RH dz DVρw Rv Ta Tr

ð13Þ

n

limits to ðDva Þ with n = 0.58 (Gat, va

2000). The molecular diffusivities are equal to 0.9691 and 0.9839 for 18O and 2H respectively (Cappa et al., 2003). The term RH refers to the relative humidity of the ambient air. Under the assumption that vapor is an ideal gas, the vapor densities and the relative humidity may be expressed as, ρva ðT a Þ ¼ eRvaðθT aa Þ, ρvr ðT r Þ ¼ eRvrðθT rr Þ ¼ v

esat ðθr Þ Rv T r

and relative humidity vertical profiles. The vertical variation of the air temperature in the atmosphere follows the average moist adiabatic lapse rate. This temporal raindrop mass change is called evaporation rate and can be expressed using the mass transfer equation,

v

eva ðθa Þ and RH ¼ ρρva ðTðTa Þa Þ ¼ eva;sat ðθa Þ , where eva(θa), evr(θr), esat(θr),

The evaporation process occurs when raindrops and the surrounding air are not in thermodynamic equilibrium. If the raindrop temperature is lower than that of the ambient air, a conductive heat flux from the air to the raindrop is detected while the raindrop looses (latent) heat since its size is reduced. The total rate of temperature change at the raindrop's surface is determined via the heat balance for the evaporating raindrop (Pruppacher and Klett, 1997),

va;sat

eva,sat(θa)are the vapor pressure and the corresponding saturated values of the ambient atmosphere and the falling raindrop at temperatures θa and θr respectively, and R∗v is the universal gas constant for the water vapor. Since the air is close to saturation at the surface of the raindrop, the vapor concentration is equal to its saturated value at the droplet temperature. The curvature effect depends on the drop size and is more pronounced for smaller drops (i.e. cloud droplets). Thus it is not included in the calculation of the saturation vapor pressure for raindrops. During their descent from the cloud base towards the ground, the raindrops fall at their terminal velocity. Due to its strong dependence on the raindrop size, the terminal velocity is higher for larger drops. Considering a constant updraft velocity w, the overall raindrop velocity is dz ¼ V ¼ V r −w dt

ð10Þ

dQ dm dh ¼ −L þ dt dt dt

ð14Þ

with L the latent heat of vaporization and dh/dt is defined in Eq. (16). The left side of Eq. (14) is defined as the total heat rate gain of the raindrop and it is responsible for the temperature difference between the raindrop and the surrounding air. The above heat gain rate can be expressed as dQ dT πρw cw VD3 dT πρw cw VD3 dT r ¼ mw cw ¼ ¼ 6 6 dz dt dt dz

ð15Þ

with cw the specific heat of water. In the right side of Eq. (14), the first term is the latent heat due to mass diffusion, while the second term refers to the sensible heat derived from the heat transfer equation (Tardif and Rasmussen, 2010),  dm 2πLDð F v Dva Þ eva;sat ðθa Þ esat ðθr Þ RH − ¼  dt Rv Ta Tr dh ¼ 2πDð F h ka ÞðT r −T a Þ: dt

with Vr the raindrop velocity and w the updraft velocity. The nomenclature for the parameters used as well as their units and mathematical expressions are given in Table A1 of the Appendix A. Based on the discussion above and Eq. (1), Eq. (9) is rewritten as,

L

dδil 3esat ðθr Þ ¼ ð F v Dva Þ dz ρw Rv D2 T r V  1 0  i !n   e δil þ 103 Dva @ i va;sat ðθa ÞT r 3 A δva þ 10 RH −  Dva esat ðθr ÞT a aieq;l−v    e va;sat ðθa ÞT r − δil þ 103 RH −1 esat ðθr ÞT a

The parameters fh and ka are the heat ventilation coefficient and the thermal conductivity of air (Table A1). Starting from Eq. (14) the differential equation that describes the vertical variation of raindrop temperature is given by,

ð11Þ

Eq. (11) describes the vertical variation of the isotopic composition of falling raindrops due to evaporation in the sub-cloud layer. The isotopic composition varies with the meteorological and isotopic conditions of the ambient atmosphere (temperature, relative humidity and isotopic composition of the atmospheric water vapor) and the raindrop characteristics (raindrop size and temperature). 3.3.1. Vertical variations of raindrop size and temperature The raindrop mass is modified due to in-flight evaporation. The raindrop diameter changes with height following the temperature

dT r 12F ka ¼ 2 h dz D ρw cw V

   F v Dva L eva;sat ðθa Þ esat ðθr Þ ðT r −T a Þ þ − RH :  F h ka Rv Ta Tr

ð16Þ

ð17Þ

3.3.2. Vertical profile of the isotopic composition of water vapor and relative humidity Schlesinger et al. (1988) studied the sub-cloud evaporation effect in precipitation using two limiting conditions, (a) constant relative humidity and (b) constant specific humidity in the sub-cloud layer. These two assumptions have important disadvantages when studying the isotopic processes occurring during precipitation. In the lower tropospheric levels, Rayleigh distillation is the dominant process that controls the water vapor isotopic composition while the mid-tropospheric vapor is mainly influenced isotopically by the non-fractionating air mass mixing

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

between vapor masses with distinct isotopic signatures (Yoshimura et al., 2011). According to the Rayleigh distillation formula δiv −δiv0 ≈

ln

Riv Riv0

!

  ¼ aeq −1 ln

  q q0

ð18Þ

the isotopic composition of water vapor depends on the initial isotopic composition of water vapor, the condensation temperature and the rainout fraction. Neglecting all other contributions that affect the isotopic variation of water vapor such as the plant transpiration, the air mass mixing, the soil evaporation etc., and according to the second assumption - constant sub-cloud specific humidity - the rainout fraction limits to unity and the isotopic composition at the ground level is the same to that of the cloud base. On the other side, the relative humidity can be altered by changing the amount of water vapor or the temperature in the atmosphere (Peixoto and Oort, 1996) explaining that a constant relative humidity (assumption (a)) is also a coarse approach for the isotopic detection of sub-cloud processes in precipitation. Under the assumption that the relative humidity is the same at all tropospheric levels, Sonntag et al. (1983) and Deshpande et al. (2010) found an exponential relationship with altitude for the isotopic composition of water vapor in the lower troposphere,      z i h  −103 δiva ðzÞ ¼ δiva zg þ 103 exp − aieq;l−v −1 H

ð19Þ

with δiva(z) and δiva(zg) the isotopic compositions of water vapor at a given vertical level z and the ground surface zg respectively and H is the atmospheric scale height. The equilibrium fractionation factor aeq,l − v is expressed at the cloud base temperature. The isotopic composition of water vapor at the cloud base is expected to reach isotopic equilibrium with precipitation and significantly differs from the isotopic content of water vapor and precipitation at the ground level (Eq. (19)). This deviation diminishes at low temperatures due to high equilibrium fractionation factors and also at higher temperatures due to secondary processes, such as below cloud evaporation and isotopic exchange. The sub-cloud evaporation effect reduces the raindrop size at each vertical level, depleting also the surrounding water vapor since the lighter isotopologues are distributed in vapor. However the isotopic contribution of the evaporated raindrops to the surrounding water vapor is negligible in comparison to the vertical composition of water vapor developed in accordance to the tropospheric specific humidity variation. Thus in Eq. (11) the isotopic composition of water vapor is expressed by Eq. (19) neglecting the small contribution from the evaporation of raindrops. 3.3.3. Initial conditions and numerical scheme The vertical profile of the isotopic composition of falling raindrops forms a system of differential equations (Eqs. (4), (11), (13), (17)) and can be defined as an initial value problem solved numerically. Here, the LSODA algorithm (Petzold, 1983) from the package deSolve (Soetaert et al., 2010) of the R language (R Core Team, 2013) is applied. Initial values for air temperature, raindrop temperature, raindrop diameter, evaporated mass fraction and isotopic composition of raindrops at the cloud base level need to be specified. The initial isotopic composition of cloud base water vapor is calculated using Eq. (19) for z = zcb while the equilibrium fractionation factor is computed at the cloud base temperature. At cloud base, vapor and liquid phases are in isotopic equilibrium and the initial isotopic composition of raindrops, δil(zcb), is expressed using the equation,   δil ðzcb Þ ¼ aieq;l−v ðT r ðzcb ÞÞ δiva ðzcb Þ−103 −103

ð20Þ

with δiva(zcb) the isotopic composition of water vapor at the cloud base and aieq,l−v(Tr(zcb)) the equilibrium fractionation factor at the raindrop temperature. The isotopic composition of water vapor at the ground

1063

level is expressed using real measurements obtained through the GNIP/ISOHIS database where a homogenous relative humidity layer is used for the simulations of the isotopic numerical model. 4. Results and discussions 4.1. δ2Η, δ18O in precipitation/water vapor and water lines for GNIP stations Stable isotopic data for water vapor/precipitation on an event basis were obtained through the GNIP/ISOHIS database for various stations located in the Northern Hemisphere. The isotopic data were selected within almost the same time period (2000–2003) except from the station in Azores where both rainfall and vapor samples were collected in 1998–1999. Summary statistics for the isotopic composition of precipitation and water vapor are shown in Table 1. δ18Ο in precipitation or water vapor (values in parentheses) ranges between −30.22‰ (−32.93‰) and 3.18‰ (−6.17‰), and δ2Η ranges between − 235.7‰ (− 225.0‰) and 22.1‰ (− 15.4‰) showing high variability. The isotopic values of water vapor reflect the reciprocal depletion in comparison to precipitation since the heavier isotopic species are preferentially distributed in the condensed phases. The observed high isotopic variability indicates that different isotopic signatures are associated with different meteorological conditions and synoptic circulation patterns since the observed isotopic composition is an integrated product of the underlying physical processes taking place during the course of the water mass from the vapor source to the sampling area and during the precipitation process. The sampling sites are divided according to their coastal proximity into two groups (A-coastal and B-continental stations). As shown in Table 1, the isotopic content for precipitation and water vapor is more depleted for group B (Ankara, Madrid and Vienna) in comparison to group A (Azores, Lisbon and Rehovot). For the continental group the average isotopic values (δ18Ο and δ2Η) of rainfall are −7.82‰ (arithmetic)/−8.2‰ (weighted) and −54.8‰ (arithmetic)/−56.2‰ (weighted) while the average values for the coastal group are −3.35‰ (−4.47‰) and −13.4‰ (−21.2‰). On the other hand, the average δ18Ο and δ2Η for water vapor are −17.98‰ and −125.2‰ for group A and −12.78‰ and −74‰ for group B. The isotopic data in precipitation and water vapor are significantly different at the 99% confidence level for the two groups according to the Wilcoxon rank-sum test (p-value bb 0.01). The isotopic depletion in continental locations is explained by the ‘continentality effect’ i.e. the progressive rainout when the air masses are moving inland from the moisture source. Continentality is controlled by the latitude of the location and the annual temperature range. A higher annual temperature range corresponds to higher seasonal temperature amplitude at stations located at similar altitudes and latitudes which in turn affects the equilibrium fractionation factors and the initial isotopic composition of an air mass prior to the precipitation event, and consequently controls the final collected isotopic composition. A measure for the deviation from the equilibrium state is also the difference between the isotopic composition of rainfall and water vapor (Δδi = δip − δiv). If the isotopic difference is comparable to the equilibrium enrichment (εi , ∗ = 103(aieq , l − v − 1)) then liquid and vapor are in isotopic equilibrium. In order to check the statistical difference between Δδi and εi,⁎ the two samples t-test is applied. In case of deuterium the ttest results are consistent for all stations. The t-statistic is strongly negative (− 3.5 b t-statistic b − 13.1) and the corresponding pvalues are lower than 0.05 indicating that Δδ2Η differs significantly from ε2,⁎. The stations of Madrid and Vienna show that Δδ18O and ε18,⁎ are significantly different (respective t-statistic = − 2.17 and − 2.48 and p-value b 0.05) while for the remaining stations the pvalue N 0.05 and the null hypothesis of the t-test cannot be rejected. Generally the t-test results show that in rainy days the rainwater and the water vapor at the ground level are not in isotopic equilibrium since the results for the deuterium imply the presence of kinetic fractionation phenomena.

1064

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

Table 1 Summary statistics of the isotopic and meteorological measurements and Meteoric Water Lines for the GNIP stations used.* Station

Azores Ankara Lisbon Madrid Rehovot Vienna Azores Ankara Lisbon Madrid Rehovot Vienna Azores Ankara Lisbon Madrid Rehovot Vienna Azores Ankara Lisbon Madrid Rehovot Vienna

δ18Ο (‰)

Type

Precipitation

Water vapor

Precipitation & Water vapor

Precipitation

Water vapor

δ2H (‰)

Min

Mean

Max

wMean⁎

Min

Mean

Max

wMean

−9.54 −30.22 −11.85 −17.03 −9.43 −24.26 −16.99 −32.93 −17.29 −24.57 −18.39 −30.27 −9.37 −30.22 −7.95 −17.03 −8.56 −15.64 −16.99 −31.67 −17.29 −24.57 −16.99 −30.27

−2.59 (±0.24) −8.82 (±0.76) −3.16 (±0.39) −7.16 (±0.67) −4.31 (±0.39) −7.48 (±0.57) −12.09 (±0.93) −18.66 (±1.29) −13.59 (±2) −16.19 (±1.18) −12.67 (±1.04) −19.09 (±2.28) −3.17 (±0.54) −8.86 (±1.54) −2.96 (±0.62) −6.58 (±0.92) −4.29 (±0.78) −7.21 (±1.5) −12.93 (±2.22) −19.8 (±3.45) −13.13 (±2.74) −16.11 (±2.26) −13.72 (±2.51) −19.09 (±3.98)

1.47 3.18 0.37 5 1.8 2.1 −6.17 −11.75 −11.33 −10.40 −9.25 −11.59 0.64 1.00 −0.62 5.00 −0.12 0.34 −6.96 −12.64 −11.33 −11.2 −10.93 −11.59

−3.38 (±0.06) −8.58 (±0.16) −4.44 (±0.09) −7.82 (±0.15) −5.58 (±0.05) −8.19 (±0.09)

−69.1 −235.7 −80.4 −129.8 −58.5 −188.4 −88.6 −225.0 −110.2 −161.2 −118.7 −210.5 −66.9 −235.7 −54.1 −129.8 −47.0 −120.0 −88.6 −225.0 −110.2 −152.4 −115.8 −210.5

−10.1 (±0.9) −60.4 (±5.1) −14.4 (±1.8) −48.6 (±4.5) −15.9 (±1.6) −55.3 (±4.2) −47.8 (±3.7) −123.9 (±8.3) −90.4 (±13.3) −115.4 (±8.4) −83.8 (±7.7) −136.2 (±16.3) −12.5 (±2.2) −61.5 (±10.7) −12.4 (±2.7) −44.3 (±6.2) −16.5 (±3.2) −52.8 (±11) −51.0 (±8.7) −139.7 (±24.3) −88.0 (±18.3) −116.2 (±16.3) −88.6 (±18.1) −135.7 (±28.3)

12.3 13.2 12.8 12.1 22.1 14 −15.4 −19.6 −70.4 −70.5 −66.3 −72.3 12.3 2 8.9 12.1 18.7 −2.4 −16.7 −91.6 −70.4 −77.1 −71.5 −72.3

−16 (±0.4) −57.4 (±1.2) −24.4 (±0.7) −53.1 (±1.1) −23.2 (±0.4) −58.2 (±0.7)

bVEL (‰)

R2

−4.12 (±0.12) −9.42 (±0.36) −3.59 (±0.18) −7.40 (±0.24) −5.44 (±0.07) −8.58 (±0.31)

Precipitation αLMWL Azores Meteoric Lines (LMWL & VEL) Ankara Lisbon Madrid Rehovot Vienna

Water vapor bLMWL (‰)

7.59 (±0.18) 9.62 (±0.7) 7.57 (±0.1) 6.98 (±0.99) 7.64 (±0.21) 9.67 (±1.05) 7.33 (±0.17) 4.21 (±1.45) 7.35 (±0.29) 18.21 (±1.75) 7.57 (±0.14) 3.83 (±1.24)

−20.1 (±1) −63.5 (±2.8) −19.3 (±1.5) −50.2 (±1.7) −23.1 (±0.5) −60.8 (±2.1)

θ (°C)

RH (%)

18.2 9.5 N/A 9.9 14.4 12.4 19.0 12.3 N/A 9.4 17.6 10.5 17.3 8 N/A 10.2 14 9.8 17.3 8 N/A 10.2 14 9.8

N/A⁎⁎ N/A N/A 81.2 N/A 76.6 N/A N/A N/A 70.7 N/A 75.2 N/A N/A N/A 82 N/A 79.2 N/A N/A N/A 82 N/A 79.2

Intersection points

R2

αT

Δα⁎⁎⁎

αVEL

0.94 0.98 0.96 0.94 0.87 0.94

8.65 8.7 N/A 8.71 8.77 8.56

1.06/0.41 1.13/0.43 N/A/0.36 1.39/0.68 1.42/0.65 0.99/0.43

4.50 (±0.27) 6.88 (±3.37) 0.62 −0.89 (±0.07) 2.9 (±0.2) 6.76 (±0.08) −4.95 (±1.44) 0.97 −14.80 (±1.02) −105.2 (±7.3) 5.65 (±0.57) −13.62 (±7.81) 0.69 −11.73 (±1.73) −79.9 (±11.8) 6.72 (±0.18) −6.69 (±2.88) 0.89 −17.95 (±1.30) −126.7 (±9.2) 5.79 (±0.26) −10.8 (±3.33) 0.81 −18.57 (±1.71) −127.3 (±9.3) 7.06 (±0.1) −1.47 (±2.01) 0.99 −10.22 (±1.22) −73.6 (±8.8)

δ18Ο (‰)

δ2Η (‰)

⁎ w-Weighted, ⁎⁎N/A- Not Available, ⁎⁎⁎Δα = αΤ − αMWL or Δα = αC − αMWL = 8 − αMWL (in the brackets).

The Meteoric Water Line for each station is determined by applying the Precipitation Weighted Least Squares Regression (Hughes and Crawford, 2012). The precipitation weighting scheme balances in some extent the contribution of the small precipitation events into the calculation of the LMWL. The ordinary least squares method is used for the determination of the Vapor Evaporation Lines. All stations have MWL slopes lower than 8 indicating the existence of secondary evaporation phenomena beneath the cloud base (Table 1). During sub-cloud evaporation more 18Ο is distributed in the liquid phase reducing the slope (Δδ2Η/Δδ18Ο) of the MWL. In order to check the statistical significance of the difference between the MWLs and the GMWL the following procedure is applied. First δ2HGMWL data are calculated based on the δ18Ο values for each station and the GMWL equation. Then the ANOVA (Analysis of Variance) and the F-ratio tests are applied taking as null hypothesis that δ2H-δ2ΗGMWL = C with C a constant while the alternative hypothesis is δ2H-δ2ΗGMWL = f(δ18Ο). For the station in Rehovot the above procedure is performed using the Eastern Mediterranean Water Line (δ2Η = 8 δ18Ο + 22‰; Gat and Carmi, 1970). The F-ratio ranges between 4.3 and 19.3 with p-values b 0.05. This reveals that the isotopic composition in each station is far from equilibrium and the GMWL or the EMWL cannot explain accurately the isotopic variability. The slopes of the MWLs are compared against the theoretical Rayleigh slopes, αT. All αLMWL values are significantly lower than αT suggesting that precipitation suffers from raindrop evaporation. Rehovot depicts the higher slope difference (aT = 8.77, Δα = α-αΤ = 1.42) indicating the highest deviation from the equilibrium state and the strongest presence of kinetic fractionation. The intercepts vary with low values in Madrid and Vienna, and high values in Rehovot while the intercepts in the remaining stations are close to the global average (10‰). In the case of Rehovot, the intercept

value (18.21‰) reflects the moisture characteristics of the Eastern Mediterranean Basin where the evaporation processes occur within low humidity continental air masses. Thus a comparison with an intercept of 10‰ may give erroneous results for the interpretation of the attributed atmospheric processes. Generally, when the intercept is higher than 10‰ the main factor that contributes to this positive difference is the continental recycling effect (Gat and Matsui, 1991); additional recycled moisture from land feeds the precipitation process altering the intercept of the MWL. Finally, the comparison of the intercept value (18.21‰) with the Eastern Mediterranean average (22‰), it suggests that secondary evaporation phenomena may contribute to the precipitation in Rehovot. The parameters of the Vapor Evaporation Lines are also shown in Table 1. The δ-pairs are placed on the left side of the GMWL (Fig. S1) exhibiting lower slopes in comparison to the LMWLs and the GMWL (Table 1). The slope of the VEL is an indicator of relative humidity conditions occurring during evaporation. The slopes range between 4.50 and 6.72 indicating evaporation under different humidity regimes (RH ≥ 50%). According to Table 1, the Azores station has the lowest slope (αVEL = 4.5) and evaporation at 50% relative humidity (Gonfiantini, 1986). However, the relatively low coefficient of determination (R2 = 0.69) indicates that the Ordinary Least Squares method is unable to describe with high predictive power the δrelationship since the vapor samples are not solely placed on a straight line but group at an elliptical scheme. The δ18Ο and δ2Η values of water vapor are distributed along a straight line in Ankara and Vienna providing high R2 values and low standard errors for the VELs parameters (Table 1). Thus the VELs for those stations can be applied to fill the missing values in the vapor isotopic time series when one of δ18Ο or δ2H has not been analyzed.

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

Further information from the Water Lines can be extracted through the intersection between the MWL and VEL. The intersection point determines the ‘initial’ isotopic composition and it is considered as the starting point prior to evaporative enrichment. A varying ‘initial’ isotopic composition is detected with the highest values in the Azores (δ18O = − 0.89‰, δ2H = 2.9‰). This station receives precipitation with an ‘early’ stage rainout where the isotopic composition is modified mainly by local meteorological factors resulting to less depleted ‘initial’ isotopic values. The ‘initial’ isotopic data for the remaining stations cannot be clustered to distinct groups according to their isotopic characteristics since low values were found in Ankara, Madrid and Rehovot and intermediate isotopic values in Lisbon and Vienna (Fig. S1; Table 1). 4.2. Isotopic enrichment of falling raindrops using an isotope-evaporation model The isotopic enrichment (Δδ) and the remaining mass fraction (fr) are calculated using the isotope-evaporation model (Eq. (5)) and the station-precipitation model (Eq. (6)) for each (θi, RHi) pair and for various raindrop sizes (1 mm ≤ D ≤ 3 mm), with i = 1,…,n and n the number of measurements. Since relative humidity data are available only for Madrid and Vienna, the calculation of the isotopic enrichment based on Stewart's model is possible only for those two sampling locations. According to Table 1 and the previous discussion in Section 4.1, Madrid and Vienna belong to the continental isotope group exhibiting similar isotopic characteristics. Therefore, the two datasets can be merged to generate a new augmented dataset for further analysis. The results from the isotopic enrichment calculations are represented in Table 2 and Fig. 1. The isotopic enrichment ranges between 0.03‰ and 20.11‰ for 18O while for 2H between 0.0‰ and 70.3‰ depending on the raindrop diameter. Higher isotopic enrichment corresponds to lower raindrop sizes since the smaller raindrops are more favorable to evaporation. The remaining mass fraction is calculated between 0% and 79% where the zero values indicate that the falling raindrops remain unaffected by sub-cloud evaporation, either because they are falling under saturated environmental conditions or because they are large enough and the evaporation process is negligible. The isotopic enrichment is linearly related to the evaporated mass fraction EF = 1 − f, for both stable water isotopes (Fig. 1 and Fig. S2a, b) according to Δδ ¼ ΔðΔδÞ ΔðEFÞ EF þ b. The slopes Δδ/EF are 0.25‰/%–0.29‰/% for δ18O and 0.9‰/%–1.0‰/% for δ2H while for the combined dataset the Δδ/EF slopes spans between 0.25‰/% and 0.28‰/% for 18O and for 2H the Δδ/EF values are around 0.9‰/%. According to the lower mass and higher vapor pressure, 1H2H16O is more favorable to evaporation than 1H18 2 O leaving the remaining raindrops less enriched in deuterium. Therefore the change in respect to the ‘initial’ isotopic composition for δ18O is expected to be higher than δil ðzcb Þ−δil ðzg Þ

δ2H. To confirm this fact the quantity Ai ð%Þ ¼ 100  ð

δil;in

Þis calcu-

lated. The numerator represents the opposite of the evaporative enrichment while the denominator δil, in is the ‘initial’ isotopic composition of liquid phase prior to evaporative enrichment. An initial guess for the starting isotopic composition (δil , in) prior to evaporative enrichment is taken by the intersection between the Meteoric Water line and the Vapor Evaporation Line for each station (Section 4.1; Table 1). Higher values of Ai correspond to 18 Ο and smaller raindrop diameter which agrees with the previous discussion. The A18O and A2H values range between 0%–57% and 0%–29% respectively. In Vienna the ‘initial’ isotopic composition is relatively enriched compared to that of Madrid (Table 1) and the corresponding Ai values are higher. Also Ai is higher for smaller raindrop diameters since the difference between the isotopic composition of the cloud base and the ground gets higher. In order to compare the above results with published works, the decrease in deuterium excess d (=δ2Η-8 δ18Ο) per evaporated fraction

1065

is calculated. Following the definition of d, the change in d per evaporat2

18

Δδ H Δδ O ed fraction is derived: Δd EF ¼ EF −8 EF , ranging from − 1.1‰/% to − 1.3‰/%. The minus sign in the Δd/EF values makes sense since the evaporation process tends to reduce d. Froehlich et al. (2008) examined the decrease in the deuterium excess in the humid Alps regions. They found a similar decrease by 1.1‰ per 1% evaporated fraction for mountain and plain stations despite of the elevation difference. Similar values for the deuterium excess change (1.1–1.2‰ decrease per 1% evaporated fraction) where found also by Kong et al. (2013) in the region of Eastern Tianshan, Central Asia. As already explained, relative humidity is an excellent indicator for the sub-cloud evaporation effect. Generally low relative humidity corresponds to high evaporation rate and consequently to high isotopic enrichment. This relationship is examined by applying a statistical model; here is the ‘heat capacity’ model, for each raindrop diameter bin. The ‘heat capacity’ model follows the equation y = A + Bx + C/x2 with x the fractional relative humidity, y the isotopic enrichment for 18 O and 2H and A, B and C the model coefficients. For each station and raindrop diameter the statistical models show high predictive performance with the coefficient of determination (R2) ranging between 81% and 91% for 18O and from 86% to 94% for 2H. Similar results are obtained for the combined data set (Table 2) with high R2 values for both isotopic species (81%–88% for 18O and 86%–92% for 2H).The model coefficients decrease for higher raindrop diameters since the isotopic enrichment is reduced with the weaker presence of the sub-cloud evaporation effect. Furthermore the isotopic enrichment for each raindrop diameter is grouped in relative humidity bins ranging from 50% to 100% with a step of 10%. The relative humidity values are rounded to the closest decade i.e. the group of RH = 90% includes isotopic enrichment data within the range 85% ≤ RH b 95%. The results are summarized using boxplots of isotopic enrichment per raindrop diameter and relative humidity bin (Fig. S2). For each case the isotopic enrichment decreases exponentially with increasing diameter. In case of RH ≥ 90%, low values of Δδ18Ο and Δδ2Η are calculated ranging from 0.22‰ - 2.69‰ and 0.9‰ - 9.8‰. In case of D = 1 mm the median values are close to 1.54‰ (RH = 90%) and 0.39‰ (RH = 100%) for Δδ18Ο, and 5.8‰ (RH = 90%) and 1.45‰ (RH = 100%) for Δδ2H. Even if the raindrops are very small, the isotopic enrichment approaches zero for almost saturated environments explaining that for RH ≥ 90% the isotopic composition does not change significantly by the sub-cloud evaporation effect.

For 50% ≤ RH ≤ 80% the isotopic enrichment follows again the exponential decrease with increasing diameter; however the isotopic composition for each raindrop diameter becomes more variable. A possible explanation for this variability could be attributed to temperature. Air temperature plays a significant role in the sub-cloud evaporation. Higher temperature values enhance evaporation, reducing the remaining raindrop mass fraction, altering also the isotopic composition in falling raindrops. In order to confirm this fact Δδ18Ο is calculated for RH = 70% and θ = 10 °C, 20 °C and 30 °C for two raindrop diameters (1.2 mm and 2.8 mm). For these calculations the initial elevation is taken equal to 1500 m. For the first raindrop diameter (D = 1.2 mm) the Δδ18O values are 7.37‰, 11.16‰ and 16.79‰ respectively. This means that for a temperature increase of 10 °C, the isotopic enrichment increases approximately by 51% and 50%. On the other hand, the larger raindrop (D = 2.8 mm) provides Δδ18Ο equal to 0.58‰, 0.89‰ and 1.30‰ with an increase of 53% and 46% per 10 °C increase. Even if the increase percentage of the isotopic enrichment is almost the same for the two raindrop diameters, the values of the isotopic enrichment are significantly different for the same temperatures. In order to explain this fact the remaining raindrop mass fraction fr is examined. The fr values for D = 1.2 mm and 2.8 mm are 0.75, 0.61, 0.40 and 0.98, 0.97, 0.95 respectively, implying that the smaller raindrop sizes are more sensitive to temperature changes. The above discussion is verified by the examination of the boxplots of the

1066

Table 2 Statistical relationship between the isotopic enrichment (Δδ) and the evaporated mass fraction (EF) and relative humidity (RH) for Madrid, Vienna and the merged data set. Statistical models

Madrid-Retiro Δδ = f(EF)

Δδ2H = f(RH)

Vienna Δδ = f(EF)

Δδ18O = f(RH)

Δδ2H = f(RH)

Merged Δδ = f(EF)

Δδ18O = f(RH)

Δδ2H = f(RH)

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

Δδ18O/ΔEF Δδ2H/ΔEF Δd/ΔΕF A (±sA) B (±sB) C (±sC) R2 A (±sA) B (±sB) C (±sC) R2

0.27 0.98 −1.14 5 (±3) −9 (±3) 3.21 (±0.45) 0.91 22 (±8) −35 (±8) 12 (±1) 0.94

0.28 1.00 −1.20 6 (±2) −8 (±2) 1.5 (±0.3) 0.89 24 (±5) −30 (±5) 5.27 (±0.92) 0.93

0.28 1.01 −1.24 5 (±1) −6 (±1) 0.8 (±0.2) 0.88 18 (±4) −22 (±3) 2.78 (±2.61) 0.92

0.28 1.01 −1.25 3.57 (±0.82) −4.17 (±0.74) 0.48 (±0.14) 0.88 14 (±2) −16 (±2) 1.67 (±0.42) 0.91

0.28 1.01 −1.27 2.65 (±0.59) −3.05 (±0.53) 0.3 (±0.1) 0.87 10 (±2) −11 (±2) 1.1 (±0.3) 0.91

0.29 1.02 −1.27 2.00 (±0.43) −2.28 (±0.39) 0.22 (±0.07) 0.87 8 (±1) −9 (±1) 0.76 (±0.22) 0.91

0.29 1.02 −1.28 1.53 (±0.33) −1.7 (±0.3) 0.16 (±0.06) 0.87 5.80 (±0.99) −6.52 (±0.89) 0.55 (±0.17) 0.91

0.29 1.02 −1.28 1.20 (±0.26) −1.36 (±0.23) 0.12 (±0.04) 0.87 4.53 (±0.76) −5.08 (±0.69) 0.41 (±0.13) 0.91

0.29 1.02 −1.28 0.9 (±0.2) −1.08 (±0.18) 0.09 (±0.03) 0.87 3.6 (±0.6) −4.02 (±0.55) 0.3 (±0.1) 0.91

0.29 1.02 −1.28 0.77 (±0.16) −0.87 (±0.15) 0.07 (±0.03) 0.87 2.90 (±0.48) −3.24 (±0.44) 0.25 (±0.08) 0.90

0.29 1.02 −1.29 0.63 (±0.13) −0.71 (±0.12) 0.06 (±0.02) 0.87 2.36 (±0.39) −2.64 (±0.36) 0.20 (±0.08) 0.90

Δδ18O/ΔEF Δδ2H/ΔEF Δd/ΔΕF A (±sA) B (±sB) C (±sC) R2 A (±sA) B (±sB) C (±sC) R2

0.25 0.87 −1.10 1 (±5) −7 (±4) 5.40 (±0.79) 0.89 1 (±14) −22 (±13) 20 (±2) 0.92

0.26 0.90 −1.19 5 (±3) −8 (±3) 2.47 (±0.55) 0.86 20 (±10) −30 (±9) 8.6 (±2) 0.90

0.27 0.91 −1.24 5 (±2) −7 (±2) 1.32 (±0.37) 0.84 19 (±7) −24 (±6) 4 (±1) 0.88

0.27 0.92 −1.26 4 (±2) −5 (±1) 0.79 (±0.26) 0.83 15 (±5) −18 (±4) 2.60 (±0.75) 0.88

0.28 0.92 −1.28 3 (±1) −4 (±1) 0.51 (±0.19) 0.83 11 (±3) −14 (±3) 1.67 (±0.54) 0.87

0.28 0.92 −1.28 2.38 (±0.83) −2.88 (±0.76) 0.36 (±0.14) 0.82 9 (±2) −10 (±2) 1.1 (±0.4) 0.87

0.28 0.93 −1.29 1.86 (±0.63) −2.23 (±0.58) 0.26 (±0.10) 0.82 7 (±2) −8 (±2) 0.8 (±0.3) 0.87

0.28 0.93 −1.29 1.47 (±0.49) −1.75 (±0.45) 0.19 (±0.08) 0.82 5 (±1) −6 (±1) 0.62 (±0.23) 0.86

0.28 0.93 −1.30 1.17 (±0.39) −1.39 (±0.36) 0.12 (±0.05) 0.82 4 (±1) −5 (±1) 0.47 (±0.19) 0.86

0.28 0.93 −1.30 0.95 (±0.31) −1.13 (±0.29) 0.09 (±0.04) 0.82 3.5 (±0.9) −4.01 (±0.82) 0.37 (±0.15) 0.86

0.28 0.93 −1.30 0.78 (±0.25) −0.92 (±0.23) 0.09 (±0.03) 0.81 2.83 (±0.73) −3.28 (±0.67) 0.30 (±0.12) 0.86

Δδ18O/ΔEF Δδ2H/ΔEF Δd/ΔΕF A (±sA) B (±sB) C (±sC) R2 A (±sA) B (±sB) C (±sC) R2

0.25 0.9 −1.1 4 (±3) −9 (±3) 4.52 (±0.52) 0.88 13 (±110) −32 (±9) 17 (±2) 0.91

0.26 0.9 −1.2 6 (±2) −9 (±2) 2.09 (±0.36) 0.85 23 (±7) −32 (±6) 7 (±1) 0.89

0.27 0.9 −1.2 5 (±1) −6 (±1) 1.13 (±0.24) 0.84 19 (±4) −24 (±4) 3.8 (±0.7) 0.88

0.27 0.9 −1.3 4 (±1) −5.0 (±0.9) 0.69 (±0.17) 0.83 14 (±3) −18 (±3) 2.29 (±0.49) 0.87

0.28 0.9 −1.3 3.04 (±0.72) −3.67 (±0.65) 0.45 (±0.12) 0.82 11 (±2) −13 (±2) 1.49 (±0.35) 0.87

0.28 0.9 −1.3 2.32 (±0.53) −2.77 (±0.48) 0.32 (±0.09) 0.82 9 (±2) −10 (±1) 1.03 (±0.26) 0.87

0.28 0.9 −1.3 1.79 (±0.40) −2.13 (±0.37) 0.23 (±0.07) 0.82 7 (±1) −8 (±1) 0.8 (±0.2) 0.86

0.28 0.9 −1.3 1.41 (±0.31) −1.66 (±0.28) 0.17 (±0.05) 0.82 5.2 (±0.9) −5.97 (±0.82) 0.56 (±0.15) 0.86

0.28 0.9 −1.3 1.12 (±0.25) −1.32 (±0.22) 0.13 (±0.04) 0.81 4.10 (±0.71) −4.74 (±0.65) 0.43 (±0.12) 0.86

0.28 0.9 −1.3 0.91 (±0.20) −1.07 (±0.18) 0.11 (±0.03) 0.81 3.32 (±0.57) −3.82 (±0.52) 0.28 (±0.08) 0.86

0.28 0.9 −1.3 0.74 (±0.16) −0.87 (±0.15) 0.09 (±0.03) 0.81 2.71 (±0.47) −3.13 (±0.42) 0.3 (±0.1) 0.86

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

Δδ18O = f(RH)

Raindrop diameter (mm) 1

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

1067

Fig. 1. Isotopic enrichment calculated from the isotope-evaporation model using the merged station dataset. a, b) linear relationship between Δδi and the evaporated mass fraction. c, d) Δδi versus relative humidity. The solid lines correspond to the ‘heat capacity’ models. The coefficients for the non-linear models are reduced for smaller raindrop sizes and as the raindrop size increases the scattered data are placed almost at a straight line.

isotopic enrichment in case of 50% ≤ RH ≤ 80% (Fig. S2). The isotopic enrichment looks highly variable for the smaller raindrop sizes providing higher interquartile ranges in the boxplots. 4.3. Isotopic enrichment using a numerical isotope evaporation model 4.3.1. Initial conditions The system of the differential equations used for the interpretation of the isotopic composition of falling raindrops and its vertical structure is described as an initial value problem. In order to apply a numerical method - here the LSODA numerical scheme - initial values for air temperature, raindrop temperature, raindrop diameter, isotopic composition of raindrops at the cloud base and evaporation rate have to be specified. However parameters such as the relative humidity (RH) and the stable isotopic composition of water vapor at the ground level must be provided in advance since they are required for solving the above differential equations. Also, they are not accounted as initial conditions for the differentiation problem since their values are assumed constant for the overall evaluation. The air temperature at ground level ranges between 5 °C to 30 °C with a step of 5 °C. Assuming that the air temperature follows the mean moist adiabatic lapse rate and using a mean cloud base elevation at about 1500 m, the cloud base temperatures are calculated and they are taken as the initial values of the air temperature for the differentiation problem. The relative humidity ranges from 40% to 90% with an increment of 10%. In order to investigate the vertical isotopic structure at almost saturated environmental conditions, a relative humidity equal to 95% is also applied. The initial raindrop temperature is calculated by the wet bulb temperature formula (Table A1) using the air temperature at the cloud base and the relative humidity. The initial raindrop diameter varies between 1 mm and 3 mm with an increment of 0.2 mm while the initial evaporated mass fraction is taken equal to zero.

As already explained initial values for the isotopic composition of raindrops at the cloud base as well as the isotopic composition of ground water vapor are necessary as initial conditions for the numerical simulations. In Section 4.2, a merged dataset is used for the evaluation of the isotope-enrichment model and the calculation of the isotopic enrichment of falling raindrops under various meteorological conditions and raindrop sizes. The average isotopic composition of ground water vapor is derived for different relative humidity values. Without splitting the isotopic composition of water vapor into relative humidity bins δ18Οv deviates from the normality case according to the Anderson– Darling test (p-value = 0.0001) under the 99% confidence level. Generally the average isotopic composition is −16.98‰ with a standard deviation of 3.13‰. Summary statistics for the isotopic composition of the ground level water vapor per relative humidity bin are represented in Table S1. δ18O in water vapor per relative humidity class ranges between − 16.39‰ and − 18.67‰ with a mean value of − 17.35‰ (s.d. 2.8‰). A first check whether the isotopic samples follow the normal distribution is performed using the Anderson-Darling test while the Kruskal–Wallis non-parametric test is applied for the examination of the (non-) distinction of the isotopic groups. According to Table S1, δ18Οv is close to normality if separate samples are assumed, providing p-values higher than 0.01 under the 99% confidence level while the null hypothesis (similar mean values) of the Kruskal-Wallis test cannot be rejected comparing the isotopic composition of water vapor for 50% ≤ RH ≤ 90% with p-value = 0.53. Thus the isotopic composition of ground level water vapor is coarsely assumed to range between −23‰ to −11% centered at −17‰. Following this coarse representation of the isotopic composition of the ground level water vapor, the initial isotopic composition of raindrops at the cloud base is derived through a two stage procedure. First the isotopic composition of the cloud base water vapor is specified using the Eq. (19) and then, under the assumption that liquid and vapor phases are in isotopic equilibrium at the

1068

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

raindrop's surface, Eq. (20) is applied for the estimation of the initial isotopic composition of raindrops. Overall the numerical isotopeevaporation model simulates the vertical isotopic structure of a falling raindrop for 2310 different initial conditions. 4.3.2. Numerical simulations Using the numerical procedure explicitly described in Section 3.3 and according to the initial meteorological and isotopic conditions described previously, the vertical variation of stable isotopes of falling raindrops is determined. Due to the large number of numerical simulations, three temperatures (5 °C, 15 °C and 30 °C) and three relative humidity (40%, 70% and 95%) cases are selected in order to discuss the impact of sub-cloud evaporation on the falling raindrops in terms of their isotopic composition. Fig. 2 visualizes the δ18O in falling raindrops while similar findings are expected for deuterium due to the strong correlation between the two isotopologues. For each pair of meteorological conditions the vertical distribution of stable isotopes of falling raindrops is simulated for various raindrop sizes. The solid lines and the confidence bands in Fig. 2 correspond to the mean and the standard deviation of the numerical simulations, since different isotopic compositions of ground level water vapor and raindrops in the cloud base are used. Generally for each single image the isotopic composition is simulated for five isotopic values of ground water vapor following the three-sigma rule described in the previous section. Then the mean and the standard deviation of the simulations are calculated for every vertical step. The isotopic enrichment is also calculated in order to investigate the contribution of sub cloud evaporation during the descent of the raindrops from the cloud base towards the ground. Table 3 provides the δ18O and δ2H enrichment in falling raindrops based on the simulations of the isotopic-numerical model. The isotopic enrichment is defined as the difference between the initial and the final isotopic composition. Also the falling distance is derived as the overall distance covered by the falling raindrops during their trip from the cloud base to the ground surface. The cloud base elevation is taken at around 1500 m above surface. Therefore falling distances lower than 1500 m imply that sub cloud evaporation processes occurring beneath the cloud base can completely evaporate the falling raindrops before reaching the surface. Generally the isotopic enrichment is higher for lower humidity and higher temperature values and for smaller raindrops. Fig. 2 shows the vertical isotopic composition of raindrops for the selected meteorological conditions. The three first plots (Fig. 2a, b and c) correspond to RH = 40% and temperatures equal to 5 °C, 15 °C and 30 °C respectively. In the case of θ = 5 °C almost all the raindrops reach the ground surface with enriched isotopic signatures. Since the sub-cloud evaporation effect is more pronounced on smaller raindrops, raindrops with D ≤ 1.2 mm are completely evaporated before reaching the ground surface and their δ18O is up to 20‰ explaining the high isotopic enrichment in falling raindrops during their descent under low humidity conditions. Since higher temperatures enhance the sub cloud evaporation effect, raindrops with D ≤ 2.2 mm cannot reach the ground surface when θ = 15 °C; when θ = 30 °C they are completely evaporated after 1000 m beneath cloud base and their isotopic composition is highly enriched. The above findings are also reported in Table 3 in terms of the isotopic enrichment and the overall traveling distance of the falling raindrops. The values in Table 3 correspond to the total isotopic enrichment of the falling raindrops and not the relative isotopic enrichment per falling distance. In case of RH = 40% and θ = 15 °C, Δδ18Ο is 32.72 ± 0.15‰ and 29.52 ± 0.13‰ for D = 2.2 mm and 1.4 mm respectively. Although Δδ18Ο is higher for D = 2.2 mm, the traveling distance of the raindrop is longer, namely 500 m more than the raindrop with D = 1.2 mm. This indicates that the kinetic fractionation is more intense for smaller raindrops; actually the induced isotopic enrichment of falling raindrops is comparable with that of the biggest raindrop (D = 2.2 mm) for a smaller falling distance. The relationship between Δδ18O and raindrop size is

more clearly depicted for those raindrops that reach the ground surface. For the same meteorological conditions as above but for D ≥ 2.4 mm, Δδ18O is reduced with raindrop size. The same behavior for the isotopic enrichment is observed for the intermediate humidity conditions of 70%. The isotopic composition of falling raindrops is significantly enriched in comparison to the initial isotopic composition (Fig. 2d, e and f). At low temperatures (θ = 5 °C) all raindrop sizes reach the ground surface whereas the isotopic variability is higher for the smaller raindrop sizes. The isotopic composition for raindrops with D ≥ 2 mm follows an almost linear relationship with falling distance and the below cloud evaporation effect does not play a significant role on the evolution of the stable isotopic composition. At higher temperatures, a stronger isotopic enrichment is observed. In the case of θ = 15 °C, raindrops with D ≤ 1.2 mm are completely evaporated covering a distance of about 1300 m below cloud base. Almost all raindrops (D ≤ 2.6 mm) do not reach the ground surface when θ = 30 °C resulting to relatively high isotopic differences compared to the isotopic composition at the cloud base. Also for the same temperature conditions but for different relative humidities, the isotopic composition increases in drier atmospheres while the falling distance decreases. Finally at RH = 95% the ambient air is almost saturated and no significant variations of the isotopic composition of falling raindrops are observed. (Table 3; Fig. 2g, h and i). Δδ18O is very low and this demonstrates that at almost saturated atmospheric conditions the sub cloud evaporation effect is negligible. The differences between the isotopic enrichment for D = 1 mm and 3 mm at the three temperature cases are 0.64‰, 0.92‰ and 1.11‰ showing that the isotopic difference is very small independently from the raindrop size's increase. The temperature dependence of the isotopic enrichment is clearly shown but Δδ18O values are comparable for each raindrop size. For example when D = 1 mm, Δδ18O is equal to 0.83‰, 1.11‰ and 1.33‰ for θ = 5 °C, 15 °C and 30 °C respectively. Even in higher temperatures the isotopic composition of the falling raindrops is hardly perturbed by the sub cloud evaporation effect and the final isotopic composition does not change significantly compared to the cloud base isotopic composition (Table 3). For the selected meteorological conditions the isotopic enrichment is also derived via the isotope-evaporation model described in Section 3.2. This model calculates the isotopic enrichment based on the prevailing meteorological conditions, without taking into account the vertical variations of temperature and the isotopic composition. Those facts are imprinted in the Δδ18Ο values leading to large differences in comparison to the numerical simulations of Δδ18O (Table 3). The nature of the isotope-evaporation model (two member mixing model) describes those significant deviations while only in the case of the near-saturated atmospheric conditions the two isotopic enrichments are generally comparable (Table 3). 5. Conclusions In this study the sub-cloud evaporation effect on precipitation is investigated in terms of its isotopic composition. In order to detect the variability of the isotopic signature due to the sub-cloud evaporation effect, three isotopic approaches are followed namely that of Meteoric Water Lines (MWLs), an isotope-evaporation model and a numerical isotope-evaporation model. First, the slope and intercept of the MWL as well as the parameters of the Vapor Evaporation Lines are examined in order to detect the presence of the bulk evaporation effect on the isotopic composition of precipitation for various stations located in the Northern Hemisphere. The deviation of the equilibrium slope is also calculated to assess the intensity of the evaporation phenomena. All stations provide slopes lower than 8 revealing the existence of kinetic fractionation beneath the cloud base whereas the lowest slopes are represented in Madrid, Spain (αLMWL = 7.33) and Rehovot, Israel (αLMWL = 7.35). Furthermore, the highest deviation from the equilibrium state is found for the Rehovot, Israel station (Δα = 1.42).

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

1069

Fig. 2. Vertical structure of the isotopic composition of falling raindrops for various meteorological conditions and raindrop sizes. a, b, c) RH = 40% T = 5 °C, 15 °C, and 30 °C. d, e, f) RH = 70% T = 5 °C, 15 °C, and 30 °C. g, h, i) RH = 95% T = 5 °C, 15 °C, and 30 °C.

In order to measure the impact of the sub-cloud evaporation effect the isotopic enrichment is calculated using an isotopic-evaporation model. This model takes into account the concurrent meteorological conditions (temperature and relative humidity) and the initial raindrop size neglecting the vertical structure of temperature and the isotopic composition of precipitation/water vapor. Due to the limited availability of meteorological information the datasets from Madrid, Spain and Vienna, Austria are merged for further analysis. Both isotopic species show a linear relationship with the evaporated mass fraction with slopes ranging between 0.25‰/% and 0.28‰/% for 18O while for 2H the

Δδ/EF values are around 0.9‰/%. Also the change of d per evaporated fraction spans −1.1‰/% and −1.3‰/% and depicts similar values with those presented in the literature (Froehlich et al. 2008; Kong et al., 2013). The isotopic enrichment is highly correlated with the raindrop diameter and relative humidity conditions. Generally, Δδ increases with decreasing relative humidity and decreasing raindrop size. Δδ is related to RH through a ‘heat capacity’ model providing high coefficients of determination (81%–88% for 18O and 86%–92% for 2H) indicating the fact that RH is an excellent indicator for the determination of the subcloud evaporation phenomenon. Temperature also contributes to the

1070 Table 3 Isotopic enrichment and falling distance for raindrops based on the numerical isotopic model. The isotopic enrichment is also provided using the isotope-mixing evaporation model. The subscripts (n) and (m) correspond to the isotopic numerical and the mixing evaporation model respectively. Cases

Raindrop diameter (mm)

RH (%)

θ (°C)

40

5

15

70

5

15

30

95

5

15

30

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

28.81 (±0.13) 129.1 (±4.6) 1150 13.77 55 22.3 (±0.1) 123 (±4.4) 650 23.47 88.7 23.60 (±0.11) 89.6 (±3.2) 290 N/A N/A 9.24 (±0.04) 64.8 (±2.3) 1500 4.7 18.6 11.13 (±0.05) 73.6 (±2.6) 1160 6.87 25.3 8.33 (±0.04) 47.9 (±1.7) 540 10.57 35.2 0.83 29 (±1) 1500 0.31 1.2 1.11 37.2 (± − 1.3) 1500 0.8 2.9 1.33 (±0.01) 38.9 (±1.4) 1500 1.37 4.5

31.97 (±0.14) 142.5 (±5.1) 1380 7.62 29.9 24.54 (±0.11) 119 (±4.2) 810 12.36 45.2 21.8 (±0.1) 82.3 (±2.9) 360 24.24 80.3 7.58 (±0.03) 52.1 (±1.9) 1500 2.86 11.2 12.04 (±0.05) 74.9 (±2.7) 1400 4.32 15.7 8.93 (±0.04) 52.6 (±1.9) 680 7.4 24 0.66 25.5 (±0.9) 1500 0.19 0.7 0.88 33.6 (±1.2) 1500 0.52 1.9 1.15 (±0.01) 37.4 (±1.3) 1500 1.12 3.6

23.3 (±0.1) 113 (±4) 1500 4.7 18.3 29.52 (±0.13) 134 (±4.8) 970 7.5 27.1 24.13 (±0.11) 92.3 (±3.3) 440 13.88 44.8 6.47 (±0.03) 44 (±1.6) 1500 1.84 7.2 10.91 (±0.05) 74 (±2.6) 1500 2.83 10.2 9.35 (±0.04) 56 (±2) 810 5.09 16.3 0.54 22.5 (±0.8) 1500 0.12 0.5 0.72 30.5 (±1.1) 1500 0.35 1.3 0.99 35.9 (±1.3) 1500 0.86 2.8

17.55 (±0.08) 82.9 (± − 2.9) 1500 3.12 12.1 28.31 (±0.13) 128 (±4.6) 1100 4.94 17.7 23.87 (±0.11) 113 (±4) 510 8.93 28.5 5.65 (±0.03) 38.2 (±1.4) 1500 1.25 4.8 9.49 (±0.04) 62.9 (±2.2) 1500 1.94 6.9 9.63 (±0.04) 58.1 (±2.1) 930 3.58 11.4 0.46 20.1 (±0.7) 1500 0.08 0.3 0.59 28 (±1) 1500 0.25 0.9 0.84 34.4 (±1.2) 1500 0.65 2.1

14.56 (±0.07) 68 (± − 2.4) 1500 2.18 8.4 29.59 (±0.13) 135.1 (±4.8) 1230 3.43 12.3 27.44 (±0.12) 107.2 (±3.8) 590 6.14 19.5 5.02 (±0.02) 33.8 (±1.2) 1500 0.89 3.4 8.47 (±0.04) 55 (±2) 1500 1.38 4.9 10.05 (±0.05) 61.4 (±2.2) 1050 2.58 8.2 0.39 18.2 (±0.7) 1500 0.06 0.2 0.50 25.8 (±0.9) 1500 0.18 0.6 0.71 32.9 (±1.2) 1500 0.49 1.6

12.54 (±0.06) 58.2 (±2.1) 1500 1.58 6.1 30.13 (±0.14) 138.2 (±4.9) 1350 2.48 8.9 28.04 (±0.13) 110.1 (±3.9) 660 4.42 14 4.52 (±0.02) 30.4 (±1.1) 1500 0.65 2.5 7.66 (±0.03) 49.8 (±1.8) 1500 1.01 3.6 10.33 (±0.05) 63.6 (±2.3) 1160 1.92 6.1 0.34 16.6 (±0.6) 1500 0.04 0.2 0.42 23.9 (±0.9) 1500 0.13 0.5 0.60 31.5 (±1.1) 1500 0.38 1.2

11.06 (±0.050) 51.1 (±1.8) 1500 1.18 4.6 32.72 (±0.15) 152.6 (±5.4) 1470 1.86 6.6 30.16 (±0.14) 120 (±4.3) 730 3.3 10.4 4.11 (±0.02) 28 (±1) 1500 0.49 1.9 7.00 (±0.03) 45.2 (±1.6) 1500 0.77 2.7 10.71 (±0.05) 66.6 (±2.4) 1270 1.46 4.6 0.3 15.2 (±0.5) 1500 0.03 0.1 0.35 22.2 (±0.8) 1500 0.1 0.4 0.50 30.2 (±1.1) 1500 0.29 0.9

9.91 (±0.04) 45.6 (±1.6) 1500 0.91 3.5 24.52 (±0.11) 109.8 (±3.9) 1500 1.43 5.1 27.51 (±0.12) 108.1 (±3.8) 790 2.52 8 3.77 (±0.02) 25.2 (±0.9) 1500 0.38 1.5 6.44 (±0.03) 41.4 (±1.5) 1500 0.59 2.1 10.95 (±0.05) 68.4 (±2.4) 1370 1.13 3.6 0.26 14.1 (±0.5) 1500 0.02 0.1 0.30 20.8 (±0.7) 1500 0.08 0.3 0.42 29 (±1) 1500 0.23 0.7

8.98 (±0.04) 41.2 (±1.5) 1500 0.72 2.8 20.62 (±0.09) 90.9 (±3.2) 1500 1.12 4 30.69 (±0.14) 122.9 (±4.4) 860 1.98 6.2 3.48 (±0.02) 23.2 (±0.8) 1500 0.3 1.1 5.97 (±0.03) 38.2 (±1.4) 1500 0.47 1.7 11.27 (±0.05) 70.9 (±2.5) 1470 0.89 2.8 0.23 13.1 (±0.5) 1500 0.02 0.1 0.26 19.5 (±0.7) 1500 0.06 0.2 0.34 28 (±1) 1500 0.18 0.6

8.22 (±0.04) 37.6 (±1.3) 1500 0.57 2.2 18.13 (±0.08) 79.1 (±2.8) 1500 0.9 3.2 29.56 (±0.13) 117.8 (±4.2) 920 1.58 5 3.23 (±0.01) 21.5 (±0.8) 1500 0.24 0.9 5.56 (±0.03) 35.5 (±1.3) 1500 0.38 1.3 10.65 (±0.05) 66.1 (±2.4) 1500 0.72 2.3 0.21 12.3 (±0.4) 1500 0.02 0.1 0.22 18.4 (±0.7) 1500 0.05 0.2 0.27 27 (±1) 1500 0.15 0.5

7.57 (±0.03) 34.6 (±1.2) 1500 0.46 1.8 16.28 (±0.07) 70.6 (±2.5) 1500 0.73 2.6 29.64 (±0.13) 118.3 (±4.2) 980 1.28 4 3.01 (±0.01) 20.1 (±0.7) 1500 0.19 0.7 5.20 (±0.02) 33.2 (±1.2) 1500 0.31 1.1 10.02 (±0.05) 61.4 (±2.2) 1500 0.59 1.8 0.19 11.5 (±0.4) 1500 0.01 0 0.19 17.4 (±0.6) 1500 0.04 0.1 0.22 25.7 (±0.9) 1500 0.12 0.4

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

30

Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰) Δδ18On (‰) Δδ2Hn (‰) Distance (m) Δδ18Om (‰) Δδ2Hm (‰)

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

1071

Table A1 Supplementary equations. Notation

Description

Equation

Unit

α18 eq,l-vO α2eq,l-vH

Equilibrium fractionation factor at temperature T (Majoube, 1971)

aeq;l−v 18 O ¼ ; expð−0:002667− 0:4156 þ 1137 Þ T T2 aeq;l−v 2 H ¼ ; expð0:052612− 76:248 þ 24844 2 Þ T T

–

Dv,a esat

Diffusivity of water vapor in air (Monteith and Unsworth,1990) Saturation vapor pressure for temperature θ (Lawrence, 2005)

P ρa

Air pressure Air density

ρw

Water density (Maidment, 1993)

þ288:9414Þðθr −3:9863Þ Þ ρw ¼ 1000ð1− ðθr508929:2ðθ r þ68:12963Þ

ka

Thermal conductivity of air (Pruppacher and Klett, 1997)

μair

ð T Þ ka ¼ 2:43  10−2 1:83210 1:71810−5 296:0

Dynamic viscosity of air (Rogers and Yau, 1989)

μ air ðTÞ ¼ 1:72 

vair

Kinematic viscosity of air

Sc, Re, Pr

Scmidt, Reynonlds and Prandtl numbers

Vr

Fall speed for raindrops (Seifert and Beheng, 2006)

Dv,a = 21.2 × 10−6(1 + 0.007T) Aθ esat ðθÞ ¼ C ; expðBþθ Þ A = 17.625, B = 243.04 °C, C = 610.94 Pa p = 101325 × (1 − 2.25577 × 10−5z)5.2559 ρa ¼ R1d T ½p−0:378RHesat ðθÞ 2

vair ¼

μ air ρa

−5 1:5 416:0 ðT−120Þ 3=2 −5 393:0 T 10 ðTþ120Þð273:15Þ

rD ; Re ¼ Vvair ; Pr ¼ vkaira Sc ¼ V r ≅aD −bD e−cD D aD = 9.65 m s−1, bD = 9.8 m s−1, 6

L = (2.501 − 0.02361θ) × 10

Latent heat of vaporization Mass ventilation coefficient (Pruppacher and Klett, 1997)

fh

Heat ventilation coefficient (Pruppacher and Klett, 1997)

1=3 1=2 Þ for ð ; Pr1=3 Re1=2 Þb1:4 f h ¼ f 1 þ 0:108ð ; Pr Re 0:708 þ 0:308ð ; Pr1=3 Re1=2 Þ for ð ; Pr1=3 Re1=2 ÞN1:4

θw

Wet bulb temperature (Stull, 2011)

θw ¼ θ ; tan−1 ½ðRH% þ 8:313659Þ1=2  þ ; tan−1 ðθ þ RH%Þ − ; tan−1 ðRH%−1:676331Þ þ 0:00391838ðRH%Þ3=2 ; tan−1 ð0:023101RH%Þ −4:686035

θd

Dew point temperature (Lawrence, 2005)

2

1=2 1=2 1=3 1=3 Þb1:4 f v ¼ f 1 þ 0:108ðSc Re Þ1=2 for ðSc Re 1=2 0:708 þ 0:308ðSc1=3 Re Þ for ðSc1=3 Re ÞN1:4 2

100

kg m−3 J m-1s−1K−1 kg m−1 s−1 m2 s−1 m s−1

cD = 600 m−1

L fv

B½ ; ln ð RH Þþ Aθ 

Pa kg m−3

–

vair Dv;a

Bþθ θd ¼ A− ; ln 100 ð RH Þ− Aθ

m2 s−1 Pa

J∙kg−1 – – °C

°C

Bþθ

A = 17.625, B = 243.04 °C

sub-cloud evaporation effect leading to higher isotopic enrichments for higher temperatures. However the temperature dependence is more distinguishable for smaller raindrops. Finally, the vertical variation of the isotopic composition of falling raindrops and the corresponding isotopic enrichment is investigated through a numerical isotope-evaporation model. This model simulates the diffusive fluxes due to evaporation in terms of the isotopic composition of the liquid phase taking into account the temperature variations of the surrounding air and the falling raindrop, the vertical change of the raindrop diameter as well as the vertical isotopic structure of the atmospheric water vapor. The isotopic composition becomes more enriched when the raindrops travel through drier and warmer atmospheres, leading to Δδ18Ο up to 20‰ depending on the raindrop size and the initial meteorological conditions. For near saturated atmospheric conditions (RH = 95%) the isotopic composition does not vary significantly indicating the absence of sub-cloud evaporation. Post condensation processes (rain evaporation and isotopic exchange) have been incorporated into the isotope enabled GCMs (Bony et al., 2008; Lee et al., 2008; Yoshimura et al., 2008; Risi et al., 2010; Dee et al., 2015) in order to describe more in detail the multiple steps of the isotope water cycle and to further improve the comparison with in-situ observations. The rain evaporation modules are mostly represented by the Stewart's model exhibiting the weaknesses of the mixing model to reproduce the ‘real’ isotopic variability. Therefore, the numerical isotope model can be an alternative of the Stewart's model in the GCMs due to its better physical representation of the sub-cloud evaporation effect and thus it could be used for validating the iso-GCMs and to perform model-station comparisons. A possible follow up of this work may be the investigation of the isotopic composition of precipitation and water vapor using high temporal resolution measurements in combination with radiosonde or satellite data for a more precise retrieval of the vertical variation of temperature and relative humidity. Together with measurements, the use of rainfall parameters such as the raindrop distribution function, the rainfall rate etc. will be helpful to derive the bulk evaporation effect in terms of the isotopic composition.

Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.scitotenv.2015.11.072.

References Blossey, P.N., Kuang, Z., Romps, D.M., 2010. Isotopic composition of water in the tropical tropopause layer in cloud-resolving simulations of an idealized tropical circulation. J. Geophys. Res. 115, D24309. Bolot, M., Legras, B., Moyer, E.J., 2013. Modelling and interpreting the isotopic composition of water vapour in convective updrafts. Atmos. Chem. Phys. 13, 7903–7935. Bony, S., Risi, C., Vimeux, F., 2008. Influence of convective proceses on the isotopic composition (δ18O and δD) of precipitation and water vapor in the tropics: 1 radiativeconvective equilibrium and tropical ocean-global atmosphere-coupled ocean–atmosphere response experiment (TOGA-COARE) simulations. J. Geophys. Res. 113, D19305. Cappa, C.D., Hendricks, M.B., DePaolo, D.J., Cohen, R.C., 2003. Isotopic fractionation of water during evaporation. J. Geophys. Res. 108 (D16), 4525. Craig, H., 1961. Isotopic variation in meteoric waters. Science 133, 1702–1703. Criss, R.E., 1999. Isotope hydrology principles of stable isotope distribution. Oxford University Press, New York. Dansgaard, W., 1964. Stable isotopes in precipitation. Tellus 16, 436–468. Dee, S., Noone, D., Buenning, N., Emile-Geay, J., Zhou, Y., 2015. SPEEDY-IER: a fast atmospheric GCM with water isotope physics. J. Geophys. Res. Atmos. 120, 73–91. Deshpande, R.D., Maurya, A.S., Kumar, B., Sarkar, A., Gupta, S.K., 2010. Rain‐vapor interaction and vapor source identification using stable isotopes from semiarid western India. J. Geophys. Res. 115, D23311. Friedman, I., Machta, I., Soller, R., 1962. Water vapour exchange between a droplet and its environment. J. Geophys. Res. 67, 2761–2766. Froehlich, K., Kralik, M., Papesch, W., Rank, D., Scheifinger, H., Stichler, W., 2008. Deuterium excess in precipitation of alpine regions-moisture recycling. Isot. Environ. Health Stud. 44 (1), 61–70. Gat, J.R., 2000. Atmospheric water balance-the isotopic perspective. Hydrol. Process. 14, 1357–2369. Gat, J., 2005. Some classical concepts of isotope hydrology. In: Aggarwal, P., Gat, J., Froehlich, K. (Eds.), Isotopes in the Water Cycle: Past, Present, and Future of a Developing Science vol. II. Springer, The Netherlands, p. 381. Gat, J.R., Carmi, I., 1970. Evolution of the isotopic composition of the atmospheric water in the Mediterranean sea area. J. Geophys. Res. 75, 3039–3048. Gat, J.R., Matsui, E., 1991. Atmospheric water balance in the Amazon Basin: an isotopic evapotranspiration model. J. Geophys. Res. 96 (D7), 13179–13188. Georgakakos, P.K., Bras, R.L., 1984. A hydrologically useful station precipitation model: 1 formulation. Water Resour. Res. 20 (11), 1585–1596. Gonfiantini, R., 1986. Environmental isotopes in lake studies. In: Fritz, P., JCh, F. (Eds.), Handbook of Environmental Isotope Geochemistry: the Terrestrial Environment B vol. II. Elsevier, Amsterdam, pp. 113–168.

1072

V. Salamalikis et al. / Science of the Total Environment 544 (2016) 1059–1072

Hughes, C.E., Crawford, J., 2012. A new precipitation weighted method for determining the meteoric water line for hydrological applications demonstrated using Australian and global GNIP data. J. Hydrol. 464-465, 42–55. IAEA/WMO, 2013. Global network of isotopes in precipitation. The GNIP database Accessible at: http://www.iaeaorg/water. Jouzel, J., 1986. Isotopes in cloud physics: multistep and multistage processes. In: Fritz, P., JCh, F. (Eds.), Handbook of Environmental Isotope Geochemistry: the Terrestrial Environment B vol II. Elsevier, Amsterdam, pp. 61–112. Klein Tank, A.M.G., et al., 2012. Daily data-set of 20th-century surface air temperature and precipitation series for the European climate assessment. Int. J. Climatol. 22, 1441–1453. Kong, Y., Pang, Z., Froenlich, K., 2013. Quantifying recycled moisture fraction in precipitation of an arid region using deuterium excess. Tellus B 65, 19251. Lawrence, M.G., 2005. The relationship between relative humidity and the dew point temperature in moist air: a simple conversion and applications. Bull. Am. Meteorol. Soc. 86, 225–233. Lee, J.E., Fung, I., 2008. ‘Amount effect’ of water isotopes and quantitative analysis of postcondensation processes. Hydrol. Process. 22 (1), 1–8. Lee, J.E., Fung, I., DePaolo, D.J., Otto, B.B., 2008. Water isotopes during the last glacial maximum: new general circulation model calaculations. J. Geophys. Res. 113, D19109. Maidment, D.R., 1993. Handbook of Hydrology. McGraw-Hill. Majoube, M., 1971. Fractionnement en 18O et en deuterium entre l'eau et sa vapour. J. Chem. Phys. 10, 1423–1428. Monteith, J.L., Unsworth, M.H., 1990. Principles of Environmental Physics. 2nd ed. Edward Arnold, London. Peixoto, J.P., Oort, A.H., 1996. The climatology of relative humidity in the atmosphere. J. Clim. 9, 3443–3463. Petzold, L., 1983. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM J. Sci. Stat. Comput. 4 (1), 136–148. Pruppacher, H.R., Klett, J.D., 1997. Microphysics of Clouds and Precipitation. 2nd ed. Kluwer Acad, Norwell Mass. R Core Team, 2013. R: A language and environment for statistical computing R Foundation for Statistical Computing Vienna Austria URL: http://wwwR-project.org/. Risi, C., Bony, S., Vimeux, F., Jouzel, J., 2010. Water-stable isotopes in the LMDZ4 general circulation model: model evaluation for present-day and past climates and

applications to climatic interpretations of tropical isotopic records. J. Geophys. Res. 115, D12118. Rogers, R.R., Yau, M.K., 1989. A Short Course in Cloud Physics. 3rd ed. Pergamon Press. Schlesinger, M.E., Oh, J.H., Rosenfeld, D., 1988. A parameterization of the evaporation of rainfall. Mon. Weather Rev. 116, 1887–1895. Seifert, A., Beheng, K.D., 2006. A two-moment cloud microphysics parameterization for mixed-phase clouds part 1: model description. Meteorog. Atmos. Phys. 92, 45–66. Sodemann, H., Masson-Delmotte, V., Schwierz, C., Vinther, B.M., Wernli, H., 2008. Interannual variability of Greenland winter precipitation sources part II: effects of north Atlantic oscillation variability on stable isotopes in precipitation. J. Geophys. Res. 113, D12111. Soetaert, K., Petzoldt, T., Setzer, W.R., 2010. Solving differential equations in R: package desolve. J. Stat. Softw. 33 (9), 1–25. Sonntag, C., Munnich, K.O., Jacob, H., Rozanski, K., 1983. Variations of deuterium and oxygen‐18 in continental precipitation and groundwater and their causes. In: Street-Perrot, A. (Ed.), Variations in the global water budget. D Reidel, Dordrecht Netherlands, pp. 107–124. Stewart, M.K., 1975. Stable isotope fractionation due to evaporation and isotopic exchanges of falling water drops: applications to atmospheric processes and evaporation of lakes. J. Geophys. Res. 80, 1133–1146. Stull, R., 2011. Wet-bulb temperature from relative humidity and air temperature. J. Appl. Meteorol. Climatol. 50, 2267–2269. Tardif, R., Rasmussen, R.M., 2010. Evaporation of non equilibrium raindrops as a fog formation mechanism. J. Atmos. Sci. 67, 345–364. Worden, J., Noone, D., Bowman, K., 2007. Importance of rain evaporation and continental convection in the tropical water cycle. Nature 445, 528–532. Yoshimura, K., Kanamitsu, M., Noone, D., Oki, T., 2008. Historical isotope simulation using reanalysis atmospheric data. J. Geophys. Res. 113, D19108. Yoshimura, K., Frankenberg, C., Lee, J., Kanamitsu, M., Worden, J., Röckmann, T., 2011. Comparison of an isotopic atmospheric general circulation model with new quasiglobal satellite measurements of water vapor isotopologues. J. Geophys. Res. 116, D19118. http://dx.doi.org/10.1029/ 2011JD016035. Zhang, X.P., Xie, Z.C., Yao, T.D., 1998. Mathematical modeling of variations on stable isotopic ratios in falling raindrops. Acta Meteorol. Sin. 12, 213–220.