Rayleigh model
Rayleigh distillation model describes the evolution of a multiple-phase system in which one phase is continuously removed to the other phase through fractional distillation. Despite its simplicity, it is a powerful framework to describe in particular the isotopic enrichment or depletion as material moves between reservoirs in an equilibrium process.
Model description
The system applying the Rayleigh equation is normally an open system from which the formed material is immediately removed. For example, the evaporation of the natural water bodies, and the formation of falling precipitation can be regarded as open systems. However, the Rayleigh equation can be also applied to other systems. One such system is a closed system (or two-phase equilibrium model; Gat, (1996)), where the material removed from one reservoir accumulates in a second reservoir in such a manner that isotopic equilibrium is maintained between the two reservoirs. An example would be the condensation of vapour to droplets in a cloud (with no falling precipitation).
The isotopic enrichment or depletion by the Rayleigh process for both open and closed systems can be mathematically established by different approaches (Fig. 1). The Rayleigh equations that are used to describe isotopic evolution of the remaining reservoir under different conditions are summarized in Table 1. The isotope ratio of the accumulated removal can be derived using mass conservation of rare isotopes between the reservoir and formed removal: [math]\displaystyle{ fR_\text{remain} + (1-f)R_\text{remove,acc} = R_0 }[/math]. In the case of condensation involving both liquid and solid phase, the accumulation is a combination of the formed removal of both phases. First, the formed liquid removal is calculated following the equation above, i.e., [math]\displaystyle{ R_\text{l,acc} = (R_0-fR_v)/(1-f) }[/math]. Then, we consider a separated Rayleigh process considering only the formation of solid phase. Following the same reasoning as above, the isotope ratio of the accumulated solid removal is calculated by [math]\displaystyle{ f_iR_\text{i,v} + (1-f_i)R_\text{i,acc} = R_\text{v0} }[/math], where [math]\displaystyle{ R_\text{v0} }[/math] and [math]\displaystyle{ R_\text{i,v} }[/math] are the isotope ratio of the remaining vapour at the start and the calculating point of the solid condensation stage, respectively, and [math]\displaystyle{ f_i }[/math] the fraction of remaining vapour within the solid condensation stage. Finally and again, using mass conservation of rare isotopes, we obtain the isotope ratio of the total accumulation of the removal by [math]\displaystyle{ (1-ff_i)R_\text{acc} = (1-f)R_\text{l,acc} + f(1-f_i)R_\text{i,acc} }[/math]. The examples of the corresponding Rayleigh profiles are shown in Fig. 2.
| Conditions | Equation (R and δ notation) | Fractionation factor used | ||
|---|---|---|---|---|
| Evaporation | open system saturated | [math]\displaystyle{ R = R_0f^{(α-1)} }[/math] | [math]\displaystyle{ δ = (1+δ_0)f^{(α-1)}-1 }[/math] | [math]\displaystyle{ α = \alpha_l = \frac{R_v}{R_l} \lt 1 }[/math] |
| open system undersaturated | [math]\displaystyle{ R = R_0f^{(α-1)} }[/math] | [math]\displaystyle{ δ = (1+δ_0)f^{(α-1)}-1 }[/math] | [math]\displaystyle{ \alpha = \frac{1}{\alpha_k \alpha_l},\quad\text{where}\quad \alpha_l = \frac{R_l}{R_v} \gt 1, \quad\text{and}\quad \alpha_k = \frac{D}{D'}(1-h)+h \gt 1, \quad\text{with}\quad h = \frac{e_v}{e_l} }[/math] | |
| closed system | [math]\displaystyle{ R = \frac{R_0}{f+\alpha(1-f)} }[/math] | [math]\displaystyle{ \delta = \frac{\delta_0-(\alpha-1)(1-f)}{f+\alpha(1-f)} }[/math] | [math]\displaystyle{ α = \frac{R_v}{R_l} \lt 1 }[/math] | |
| Condensation | saturation over ice | [math]\displaystyle{ R = R_0f^{(α-1)} }[/math] | [math]\displaystyle{ δ = (1+δ_0)f^{(α-1)}-1 }[/math] | [math]\displaystyle{ α = \alpha_s = \frac{R_s}{R_v} \gt 1 }[/math] |
| supersaturation over ice | [math]\displaystyle{ R = R_0f^{(α-1)} }[/math] | [math]\displaystyle{ δ = (1+δ_0)f^{(α-1)}-1 }[/math] | [math]\displaystyle{ \alpha = \alpha_k \alpha_s, \quad\text{where}\quad \alpha_s = \frac{R_s}{R_v} \gt 1, \quad\text{and}\quad \alpha_k = \frac{S_i}{\alpha_s \frac{D}{D'}(S_i-1)+1} \lt 1, \quad\text{with}\quad S_i = \frac{e_v}{e_i} }[/math] | |
Figure 1 Schematic presentation of the mathematical approaches of the Rayleigh model in open system (a, b) and closed system (c). Adapted from Gat et al (2001). (a) From a reservoir containing N abundant isotopologues (e.g. H216O) and Ni rare isotopologues (e.g. H218O or HD16O) small amounts of both species, dN and dNi respectively, are being removed under equilibrium fractionation conditions: [math]\displaystyle{ dN_i/dN = αR }[/math]. (b) The changing isotopic composition of the reservoir is calculated from a mass balance consideration for the rare isotopologues: [math]\displaystyle{ RN = (R-dR)(N-dN) + αRdN }[/math]. (c) In closed system, the changing isotopic composition of the reservoir can be calculated from a mass balance consideration for the rare isotopologues: [math]\displaystyle{ R_0N_0 = R(N_0-dN) + αRdN }[/math], with the remaining fraction [math]\displaystyle{ f = (N_0-dN)/N_0 }[/math].
Figure 2 Isotopic change during (a), (c) Rayleigh evaporation, and (b), (d) Rayleigh condensation. Rayleigh evaporation occurs at 20 °C (thus a constant fractionation factor) for initial liquid compositions of δ18O = -10 ‰ and δD = -70 ‰ (d = 10 ‰), in open system with unsaturated (RH = 75 %; black) or saturated (grey) environment, and closed system (blue). Rayleigh condensation occurs under continuous cooling (thus also a contentiously changing fractionation factor) from T = 20 °C and RH = 75 %. The initial vapour compositions are δ18O = -13 ‰ and δD = -94 ‰ (d = 10 ‰). For the condensation to ice below 0 °C, two circumstances are presented. The saturation circumstance (grey) is a classical Rayleigh process where vapour forms ice crystals under equilibrium conditions, using the saturation pressure over ice (ei). The supersaturation circumstance (black) takes into account the supersaturation over ice where the ambient vapour pressure is ev = Si*ei. Si is the defined saturation ratio, as ev/ei, and here takes the form Si = 1-0.004T after Risi et al (2010a). In this circumstance, the fractionation factor combining equilibrium and kinetic effects given by Jouzel & Merlivat (1984) is used.
Model codes in MATLAB
Rayleigh evaporation
% rayleigh_evaporation.m
%--------------------------------
% Description:
% simple Rayleigh evaporation model, after Stern and Blisniuk, 2002
%--------------------------------
% input variables: T0 (Kelvin), d18O0 (permil), d2H0 (permil),
% system (open system = 1, closed system = 0), varargin (relative humidity)
% run example: rayleigh_evaporation(20+273.15,-10,-70,1,0.75)
% rayleigh_evaporation(20+273.15,-10,-70,1)
% rayleigh_evaporation(20+273.15,-10,-70,0)
%--------------------------------
% HS, Do 12 Dez 2013 15:43:28 CET
% YW, structured and included kinetic effect, 29 Aug 2020
%--------------------------------
function data = rayleigh_evaporation(T0,d18O0,d2H0,system,h)
% %% dubugging
% T0 = 20+273.15;
% d18O0 = -10; % close to Bergen precipitation average
% d2H0 = -70;
% system = 1;
% h = 0.75;
%% fractionation factors
%--------------------------------
% Majoube, M., 1971a: Fractionnement en oxyg`ene-18 entre la glace et la vapeur d’eau. J.
% Chim. Phys., 68, 625–636.
% Majoube, M., 1971b: Fractionnement en oxyg`ene 18 et deut ́erium entre l’eau et sa vapeur. J. Chim.
% Phys., 68, 1423–1436.
% Merlivat, L. and G. Nief, 1967: Fractionnement isotopique lors des changements d‘ ́etat
% solide-vapeur et liquide-vapeur de l’eau `a des temp ́eratures inf ́erieures a ` 0 ◦ C. Tellus,
% 19, 122–127.
%--------------------------------
function d = a2H(T) % T in K
d = 1/exp(24.844e3/(T^2)-76.248/T+52.612e-3); % Majoube 1971b, R_vapour/R_liquid < 1
end
function d = a18O(T) % T in K
d = 1/exp(1137/(T^2)-0.4156/T-0.0020667); % Majoube 1971b, R_vapour/R_liquid < 1
end
%% initialization
% isotope ratios of VSMOW reference water
R2std = 155.76e-6;
R18std = 2005.20e-6;
% initial isotope ratio
R18O0 = (d18O0/1000+1)*R18std;
R2H0 = (d2H0/1000+1)*R2std;
% % ouput variables
% data.T = [];
% data.f = [];
% data.d18Ov = [];
% data.d18Ol = [];
% data.d2Hv = [];
% data.d2Hl = [];
% data.dxsv = [];
% data.dxsl = [];
% data.d18Ov_acc = [];
% data.d2Hv_acc = [];
% data.dxsv_acc = [];
%% Rayleigh evaporation process: (1) in an open system (material is immediately removed)
% e.g., condensation of vapour to droplets and falls immediately from a cloud (precipitation occurs)
if system == 1
fprintf('Rayleigh evaporation process: in an open system (material is immediately removed) \n');
fprintf('T0 = %6.2f K (%4.2f °C) \n', T0, T0-273.15);
fprintf(' f d18Ol d18Ov d2Hl d2Hv dxsl dxsv d18Ov_acc dxsv_acc \n');
n = 1; % iteration steps
for f = 1:-0.02:0 % remaining liquid fraction
% relative humidity
if ~exist('h','var')
h = 1;
end
% modified Rayleigh evaporation model including kinectic effct
% during the vapour formation at undersaturated conditions (Jouzel and Merlivat 1984)
% Diffusivities D adopted from Merlivat 1978:
% D_H18O/D_H16O = 0.9723 ± 0.0007
% D_HD16O/D_H216O = 0.9755 ± 0.0009
ak_18O = 1/(1/0.9723*(1-h)+h);
ak_2H = 1/(1/0.9755*(1-h)+h);
aeff_18O = a18O(T0) * ak_18O;
aeff_2H = a2H(T0) * ak_2H;
% isotope ratio of remaining liquid
R18Ol = R18O0*f^(aeff_18O-1);
R2Hl = R2H0*f^(aeff_2H-1);
% (1) delta value of remaining liquid
d18Ol = (R18Ol/R18std-1)*1000;
d2Hl = (R2Hl/R2std-1)*1000;
dxsl = d2Hl-8*d18Ol;
% (2) delta value of the instantaneously removed vapour (in equilibrium with the current remaining liquid)
d18Ov = ((d18Ol/1000+1)*aeff_18O-1)*1000;
d2Hv = ((d2Hl/1000+1)*aeff_2H-1)*1000;
dxsv = d2Hv-8*d18Ov;
% (3) isotope ratio of the accumulated removed vapour
R18Ov_acc = R18O0*(1-f^aeff_18O)/(1-f); % Mook 2001, Page 38
R2Hv_acc = R2H0*(1-f^aeff_2H)/(1-f);
% delta value of the accumulated removed vapour
d18Ov_acc = (R18Ov_acc/R18std-1)*1000;
d2Hv_acc = (R2Hv_acc/R2std-1)*1000;
dxsv_acc = d2Hv_acc-8*d18Ov_acc;
% % (3) Calculating the isotope ratio of the accumulated removal using mass conservation approach
% % mass conservation of rare isotopes between the reservoir and formed removal: f*Rl+(1-f)*Rv_acc=R0,
% % after rearranging, we have: Rv_acc=(R0-f*Rl)/(1-f)
% R18Ov_acc2 = (R18O0-f*R18Ol)/(1-f);
% R2Hv_acc2 = (R2H0-f*R2Hl)/(1-f);
% data.d18Ov_acc2(n) = (R18Ov_acc2/R18std-1)*1000;
% data.d2Hv_acc2(n) = (R2Hv_acc2/R2std-1)*1000;
% data.dxsv_acc2(n) = data.d2Hv_acc2(n)-8*data.d18Ov_acc2(n);
fprintf('%8.4f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n',f,d18Ol,d18Ov,d2Hl,d2Hv,dxsl,dxsv,d18Ov_acc,dxsv_acc);
% assign to output
data.T(n) = T0;
data.f(n) = f;
data.d18Ol(n) = d18Ol;
data.d18Ov(n) = d18Ov;
data.d2Hl(n) = d2Hl;
data.d2Hv(n) = d2Hv;
data.dxsl(n) = dxsl;
data.dxsv(n) = dxsv;
data.d18Ov_acc(n) = d18Ov_acc;
data.d2Hv_acc(n) = d2Hv_acc;
data.dxsv_acc(n) = dxsv_acc;
% iterate to next step
n = n+1;
end
end
%% Rayleigh evaporation process: (2) in closed system (material maintained under two-phase equilibrium)
% i.e., the material removed stays in the closed system and stays in equilibrium with the reservoir
% e.g., condensation of vapour to droplets in a cloud (only cloud is formed, but precipitation does not occur)
if system == 0
fprintf('Rayleigh evaporation process: in an closed system (material maintained under two-phase equilibrium) \n');
fprintf('T0 = %6.2f K (%4.2f °C) \n', T0, T0-273.15);
fprintf(' f d18Ol d18Ov d2Hl d2Hv dxsl dxsv d18Ov_acc dxsv_acc \n');
n = 1; % iteration steps
for f = 1:-0.02:0 % remaining liquid fraction
% isotope ratio of remaining liquid
R18Ol = R18O0/(a18O(T0) - f*(a18O(T0)-1)); % Gat1996
R2Hl = R2H0/(a2H(T0) - f*(a2H(T0)-1)); % Gat1996
% (1) delta value of remaining liquid
d18Ol = (R18Ol/R18std-1)*1000;
d2Hl = (R2Hl/R2std-1)*1000;
dxsl = d2Hl-8*d18Ol;
% (2) delta value of the instantaneously removed vapour (in equilibrium with the current remaining liquid)
d18Ov = ((d18Ol/1000+1)*a18O(T0)-1)*1000;
d2Hv = ((d2Hl/1000+1)*a2H(T0)-1)*1000;
dxsv = d2Hv-8*d18Ov;
% (3) delta value of the accumulated removed vapour
d18Ov_acc = d18Ov;
d2Hv_acc = d2Hv;
dxsv_acc = dxsv;
fprintf('%8.4f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n',f,d18Ol,d18Ov,d2Hl,d2Hv,dxsl,dxsv,d18Ov_acc,dxsv_acc);
% assign to output
data.T(n) = T0;
data.f(n) = f;
data.d18Ol(n) = d18Ol;
data.d18Ov(n) = d18Ov;
data.d2Hl(n) = d2Hl;
data.d2Hv(n) = d2Hv;
data.dxsl(n) = dxsl;
data.dxsv(n) = dxsv;
data.d18Ov_acc(n) = d18Ov_acc;
data.d2Hv_acc(n) = d2Hv_acc;
data.dxsv_acc(n) = dxsv_acc;
% iterate to next step
n = n+1;
end
end
end
% fin
Rayleigh condensation
% rayleigh_condensation.m
%--------------------------------
% Description:
% - simple Rayleigh condensation model, after Stern and Blisniuk, 2002
% - temperature dependent fractionation factor varies during cooling process
% (1) isobaric cooling: e.g., air parcel advects to cold regions
% (2) or pseudoadiabatic cooling: e.g., updraft or lift of surface air parcel
%--------------------------------
% input variables: T0 (Kelvin), RH0 (%), d18O0 (permil), d2H0 (permil),
% sice: saturation conditions over ice when temperature <0 degree
% > 1: using the constant supersaturation given as sice;
% = 1: meaning alpha_k = 1, thus no supersaturation;
% = 0: using formula Si = 1-0.004*T (in Celsius)
% run example: rayleigh_condensation(20+273.15,0.75,-13,-94,0)
% rayleigh_condensation(20+273.15,0.75,-13,-94,1)
% rayleigh_condensation(20+273.15,0.75,-13,-94,1.2)
%--------------------------------
% HS, Do 12 Dez 2013 15:43:28 CET
% YW, structured and included kinetic effect, 29 Aug 2020
%--------------------------------
% YW, 26 Mar 2020
% limitation:
% Rayleigh condensation model only reflects the process of water vapor removing
% and ignores: mixing of air masses, cloud water-vapor isotope exchange,
% subcloud droplet evaporation, and radiative cooling or heating.
%--------------------------------
function data = rayleigh_condensation(T0,RH0,d18O0,d2H0,sice)
% %% for debugging
% T0 = 20+273.15;
% RH0 = 0.75; % Craig-Gordon, 1965 (Fig. 15)
% d18O0 = -13;
% d2H0 = -94;
% sice = 0;
% RH0 = 0.91; % 20161207 AR event at Bergen
% d18O0 = -12.5;
% d2H0 = -100;
%% functions adopted
% fractionation factors
%--------------------------------
% Majoube, M., 1971a: Fractionnement en oxyg`ene-18 entre la glace et la vapeur d’eau. J.
% Chim. Phys., 68, 625–636.
% Majoube, M., 1971b: Fractionnement en oxyg`ene 18 et deut ́erium entre l’eau et sa vapeur. J. Chim.
% Phys., 68, 1423–1436.
% Merlivat, L. and G. Nief, 1967: Fractionnement isotopique lors des changements d‘ ́etat
% solide-vapeur et liquide-vapeur de l’eau `a des temp ́eratures inf ́erieures a ` 0 °C. Tellus,
% 19, 122–127.
%--------------------------------
function d = a2H(T) % T in Kelvin
if T >= 273.15 % above 0 celsius
d = exp(24.844e3/(T^2)-76.248/T+52.612e-3); % Majoube 1971b, R_liquid/R_vapour > 1
else
d = exp(16289./T.^2 - 0.0945); % Merlivat & Nief 1967, R_solid/R_vapour > 1
end
end
function d = a18O(T) % T in Kelvin
if T >= 273.15 % above 0 celsius
d = exp(1137/(T^2)-0.4156/T-0.0020667); % Majoube 1971b, R_liquid/R_vapour > 1
else
d = exp(11.839./T - 0.028224); % Majoube 1971a, R_solid/R_vapour > 1
end
end
% supersaturation with respect to ice
% Si = 1 − λT, with λ = 0.004. From Graf 2017 dissertation, page 53, Fig 3.1
% (as in Risi et al., 2010b, after Jouzel and Merlivat 1984.)
% Risi, C., S. Bony, F. Vimeux, and J. Jouzel, 2010b: 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.
% sice: > 1: using the constant supersaturation given here;
% = 1: meaning alpha_k = 1, thus no supersaturation;
% = 0: using formula Si = 1-0.004*T (in Celsius)
function s = si(T)
if sice >= 1
s = sice;
else % elseif sice == 0
s = 1-0.004*(T-273.15);
end
end
% vapour saturation pressure
function es = esat(T) % T in Kelvin
if T >= 273.15
% es = 611*exp((17.27*T)/(273+T)); % T in celsius
es = 611*exp(17.27*(T-273.15)/T); % T in Kelvin
else
% Saturation vapour pressure over ice (Murphy and Koop, 2005) for T > 110 K:
es_ice = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T);
% the actual saturation vapour pressure
es = es_ice * si(T);
end
end
% saturated vapour mixing ratio
function ws = wsat(es,p)
ws = 0.622*es/(p-es);
end
% pressure at height of z
function p = pressure(z)
p0 = 101325; % Pa
z0 = 8000; % m scale height
p = p0*exp(-z/z0);
end
% moist adiabatic lapse rate [Curry and Webster,1999] (cited after Stern and Blisniuk, 2002)
function lr = gamma(ws,T) % T in Kelvin
% physical constants
g = 9.80665; % m/s^2 the gravitational acceleration
Cpd = 1.00464e3; % J/kg/K the specific heat capacity of dry air at constant pressure
L = 2.5104e6; % J/kg the latent heat of vaporization or sublimation
Rd = 287.04; % J/kg/K the specific gas constant for dry air
lr = g/Cpd*(1+L*ws/Rd/T)/(1+0.622*L*L*ws/Rd/Cpd/T/T); % [°C/m]
end
%% initialization
% isotope ratios of VSMOW reference water
R2std = 155.76e-6;
R18std = 2005.20e-6;
% initial isotope ratio
R18O0 = (d18O0/1000+1)*R18std;
R2H0 = (d2H0/1000+1)*R2std;
dxs0 = R2H0-8*R18O0;
% initial conditions
f = 1.0; % remaining vapour fractionation
z = 0.0; % elevation
p = pressure(z);
w0 = RH0*wsat(esat(T0),p); % kg/kg
% convert water mixing ratio unit
md = 28.97; % [g/mol] dry air molar mass
mv = 18.02; % [g/mol] water vapour molar mass
wconc0 = w0*md/mv*1e6; % [kg/kg] --> [ppmv]
rh = RH0;
dt = 0.5; % temperature step
T = T0+dt;
%% for the circumstance of isobaric cooling
%% (1) cool to saturation
fprintf('Rayleigh condensation process: in an open system (material is immediately removed) \n');
fprintf('T0 = %6.2f K (%4.2f °C) \n', T0, T0-273.15);
fprintf(' T (°C) d18Ov d2Hv dxsv rh f WConc[ppmv] \n');
n = 1; % iteration number
while rh < 0.99
T = T-dt; % decrease temperature by dt °C at each iteration
ws = wsat(esat(T),p); % saturated water mixing ratio at each iteration
rh = w0/ws; % RH at each iteration
wconc = wconc0; % [kg/kg] --> [ppmv]
fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',T-273.15,d18O0,d2H0,dxs0,rh,f,wconc);
% assign to output
data.T(n) = T;
data.f(n) = f;
data.rh(n) = rh;
data.Wconc(n) = wconc;
data.d18O(n) = d18O0;
data.dD(n) = d2H0;
data.dxs(n) = d2H0-8*d18O0;
data.d18Oc(n) = NaN; % condensate
data.d2Hc(n) = NaN;
data.dxsc(n) = NaN;
data.d18O_acc(n) = NaN; % accumulated condensate (liquid+solid)
data.d2H_acc(n) = NaN;
data.dxs_acc(n) = NaN;
% iterate to next step
n = n+1;
end
fprintf('Reached saturation now\n');
ws0 = w0;
%% (2) continue to cool and form liquid condensate
while T > 273.15 % isobaric cooling with condensation to 0 °C
T = T-dt; % decrease temperature by dt °C at each iteration
ws = wsat(esat(T),p); % saturated water mixing ratio at each iteration
wconc = ws*md/mv*1e6; % [kg/kg] --> [ppmv]
f = ws/ws0;
% isotope ratio of remaining vapour
R18Ov = R18O0*f^(a18O(T)-1);
R2Hv = R2H0*f^(a2H(T)-1);
% (1) delta value of remaining vapour
d18Ov = (R18Ov/R18std-1)*1000;
d2Hv = (R2Hv/R2std-1)*1000;
dxsv = d2Hv - 8*d18Ov;
% (2) delta value of the instantaneously removed condensate (in equilibrium with the current remaining vapour)
d18Oc = ((d18Ov/1000+1)*a18O(T)-1)*1000;
d2Hc = ((d2Hv/1000+1)*a2H(T)-1)*1000;
dxsc = d2Hc-8*d18Oc;
% (3) Calculating the isotope ratio of the accumulated removal using mass conservation approach
% mass conservation of rare isotopes between the reservoir and formed removal: f*Rv + (1-f)*Rl_acc = R0,
% after rearranging, we have: Rl_acc = (R0-f*Rv)/(1-f)
R18Ol_acc = (R18O0-f*R18Ov)/(1-f);
R2Hl_acc = (R2H0-f*R2Hv)/(1-f);
d18Ol_acc = (R18Ol_acc/R18std-1)*1000;
d2Hl_acc = (R2Hl_acc/R2std-1)*1000;
dxsl_acc = d2Hl_acc-8*d18Ol_acc;
fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',T-273.15,d18Ov,d2Hv,dxsv,rh,f,wconc);
% assign to output
data.T(n) = T;
data.f(n) = f;
data.rh(n) = rh;
data.Wconc(n) = wconc;
data.d18O(n) = d18Ov;
data.dD(n) = d2Hv;
data.dxs(n) = dxsv;
data.rh(n) = rh;
data.d18Oc(n) = d18Oc;
data.d2Hc(n) = d2Hc;
data.dxsc(n) = dxsc;
data.d18O_acc(n) = d18Ol_acc;
data.d2H_acc(n) = d2Hl_acc;
data.dxs_acc(n) = dxsl_acc;
% iterate to next step
n = n+1;
end
fprintf('Have cooled to 0 celcius now \n');
ws0_ice = ws;
R18O0 = R18Ov;
R2H0 = R2Hv;
f0 = f;
%% (3) continue to cool and form solid condensate
% considering mixed-phase cloud: subsaturated over liquid droplets and supersaturated over ice crystals
% closure assumption: liquid droplets --> evaporate to vapour --> deposite on ice crystals
% therefore: when T<0 celsius, the reduced fraction of vapour is due to its deposit on ice crystals with a fractionation
% factor alpha_eff = alpha_eq*alpha_k taking into account the kinetic effct under the supersaturation condition.
% i.e.: we start a new Rayleigh condensation model starting with the remaining vapour
% after the liquid condensation stage above, and use the new factionation factor alpha_eff
while T > -60+273.15
T = T-dt;
ws = wsat(esat(T),p);
wconc = ws*md/mv*1e6; % [kg/kg] --> [ppmv]
f = ws/ws0; % the total fraction of remaining vapour through all stages
fi = ws/ws0_ice; % the fraction of remaing vapour in the ice condensation stage only
% modified Rayleigh condensation model including kinectic effct
% during the ice formation at supersaturated conditions (Jouzel and Merlivat 1984)
% Diffusivities D adopted from Merlivat 1978:
% D_H18O/D_H16O = 0.9723 ± 0.0007
% D_HD16O/D_H216O = 0.9755 ± 0.0009
ak_18O = si(T)/(a18O(T)*1/0.9723*(si(T)-1) + 1);
ak_2H = si(T)/(a2H(T)*1/0.9755*(si(T)-1) + 1);
aeff_18O = a18O(T) * ak_18O;
aeff_2H = a2H(T) * ak_2H;
% isotope ratio of remaining vapour
R18Ov = R18O0*fi^(aeff_18O-1); % Jouzel and Merlivat, 1984
R2Hv = R2H0*fi^(aeff_2H-1); % Jouzel and Merlivat, 1984
% (1) delta value of remaining vapour
d18Ov = (R18Ov/R18std-1)*1000;
d2Hv = (R2Hv/R2std-1)*1000;
dxsv = d2Hv - 8*d18Ov;
% (2) delta value of the instantaneously removed condensate (in equilibrium with the current remaining vapour)
d18Oc = ((d18Ov/1000+1)*aeff_18O-1)*1000;
d2Hc = ((d2Hv/1000+1)*aeff_2H-1)*1000;
dxsc = d2Hc-8*d18Oc;
% (3) Calculating the isotope ratio of the accumulated removal using mass conservation approach
% (i) mass conservation of rare isotopes between the reservoir and formed removal:
% fi*Rv + (1-fi)*Ri_acc = Rv0,
% where fi is the fraction of remaing vapour in the ice condensation stage, and
% Rv0 is the isotope ratio of the remaining vapour after liquid condensation stage
% after rearranging, we have: Ri_acc = (Rv0-fi*Rv)/(1-fi)
R18Oi_acc = (R18O0-fi*R18Ov)/(1-fi);
R2Hi_acc = (R2H0-fi*R2H0)/(1-fi);
% d18Oi_acc = (R18Oi_acc/R18std-1)*1000;
% d2Hi_acc = (R2Hi_acc/R2std-1)*1000;
% dxsi_acc = d2Hi_acc-8*d18Oi_acc;
% (ii) the total accumulation is sum of the liquid (1-f0) and solid condensate (f0*(1-fi))
% (1-f0*fi)*R_acc = (1-f0)*Rl_acc + (f0*(1-fi))*Ri_acc
R18O_acc = ((1-f0)*R18Ol_acc + (f0*(1-fi))*R18Oi_acc)/(1-f0*fi);
R2H_acc = ((1-f0)*R2Hl_acc + (f0*(1-fi))*R2Hi_acc)/(1-f0*fi);
d18O_acc = (R18O_acc/R18std-1)*1000;
d2H_acc = (R2H_acc/R2std-1)*1000;
dxs_acc = d2H_acc-8*d18O_acc;
fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',T-273.15,d18Ov,d2Hv,dxsv,rh,f,wconc);
% assign to output
data.T(n) = T;
data.f(n) = f;
data.rh(n) = rh;
data.Wconc(n) = wconc;
data.d18O(n) = d18Ov;
data.dD(n) = d2Hv;
data.dxs(n) = dxsv;
data.rh(n) = rh;
data.d18Oc(n) = d18Oc;
data.d2Hc(n) = d2Hc;
data.dxsc(n) = dxsc;
data.d18O_acc(n) = d18O_acc;
data.d2H_acc(n) = d2H_acc;
data.dxs_acc(n) = dxs_acc;
% iterate to next step
n = n+1;
end
%% for the circumstance of pseudoadiabatic cooling
% - establish the relationship between T and z
z = 0;
p = pressure(z);
n = 1;
for T = T0:-dt:-60+273.15
data.T(n) = T;
data.z(n) = z;
data.p(n) = p;
% the elevation change
ws = wsat(esat(T),p);
lr = gamma(ws,T); % lapse rate
dz = dt/lr; % dt = lr*dz
% the new elevation and pressure
z = z+dz;
p = pressure(z);
data.lr(n) = lr;
% iterate to next step
n = n+1;
end
end
% fin
Example of usage
1. Generate the desired data sets
%% create data sets % Rayleigh evaporation T0 = 20+273.15; d18O0 = -10; % close to Bergen precipitation average d2H0 = -70; h = 0.75; % relative humidity de_ok = rayleigh_evaporation(T0,d18O0,d2H0,1,h); % in an open system, unsaturated (including kinetic effect) de_oe = rayleigh_evaporation(T0,d18O0,d2H0,1); % in an open system, saturated (equilibrium) de_cl = rayleigh_evaporation(T0,d18O0,d2H0,0); % in a closed system % Rayleigh condensation T0 = 20+273.15; RH0 = 0.75; % Craig-Gordon, 1965 (Fig. 15) d18O0 = -13; d2H0 = -94; % > 1: using the constant supersaturation given as sice; % = 1: % = 0: dc_ln = rayleigh_condensation(T0,RH0,d18O0,d2H0,0); % using linear formula Si = 1-0.004*T (in Celsius) dc_no = rayleigh_condensation(T0,RH0,d18O0,d2H0,1); % meaning alpha_k = 1, thus no supersaturation, e = e_ice dc_11 = rayleigh_condensation(T0,RH0,d18O0,d2H0,1.1); % using the constant supersaturation as given
2. Plot Rayleigh evaporation and condensation profiles.
As an example, Fig. 2 is created using the data sets generated from the above step and using the plotting codes from below.
%% plot Rayleigh evaporation and condensation model
fs = 14;
clear h
figure
set(gcf,'color','white','position',[0 0 1200 1150],'PaperPositionMode','auto');
% left_color = [0 0 0]; % black
% right_color = [0 0.4470 0.7410]; % blueish % blue [0 0 1];
% set(gcf,'defaultAxesColorOrder',[left_color; right_color]);
% Rayleigh evaporation
h(1) = subplot(2,2,1);
% set(h(1),'Position', [0.1 0.56 0.38 0.38]);
hold on
pe(1) = plot(de_ok.f,de_ok.d18Ol,'k-','linewidth',2);
pe(2) = plot(de_ok.f,de_ok.d18Ov,'k--','linewidth',2);
pe(3) = plot(de_ok.f,de_ok.d18Ov_acc,'k:','linewidth',2);
pe(4) = plot(de_oe.f,de_oe.d18Ol,'-','color',[0.6 0.6 0.6],'linewidth',2); %'color',[0.3010 0.7450 0.9330],
pe(5) = plot(de_oe.f,de_oe.d18Ov,'--','color',[0.6 0.6 0.6],'linewidth',2);
pe(6) = plot(de_oe.f,de_oe.d18Ov_acc,':','color',[0.6 0.6 0.6],'linewidth',2);
pe(7) = plot(de_cl.f,de_cl.d18Ol,'-','color',[0.3010 0.7450 0.9330],'linewidth',2);
pe(8) = plot(de_cl.f,de_cl.d18Ov,'--','color',[0.3010 0.7450 0.9330],'linewidth',2);
pe(9) = plot(de_cl.f,de_cl.d18Ov_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2);
ylim(h(1),[-27 40]);
xlabel('Residual water fraction','fontsize',fs);
ylabel(['\delta^{18}O / ',char(8240)],'fontsize',fs);
lgd1 = legend(h(1),pe(1:3),'residual water','instant vapour','total vapour removed','Location','northwest');
lgd1_pos = lgd1.Position;
lgd1.FontSize = fs;
% 2nd axes
a2 = axes('position',get(h(1),'position'),'visible','off'); % create the same axes but set it invisible
lgd2 = legend(a2,pe(4:6),'','','','location','northwest','box','off');
lgd2_pos = lgd1_pos+[-lgd1_pos(3)*0.304 -lgd1_pos(4)*0.72 0 lgd1_pos(4)*0.8]; % move the 2nd legend to the desired position
set(lgd2,'box','off','Position',lgd2_pos);
% 3rd axes
a3 = axes('position',get(gca,'position'),'visible','off'); % create the same axes but set it invisible
lgd3 = legend(a3,pe(7:9),'','','','location','northwest','box','off');
lgd3_pos = lgd1_pos+[-lgd1_pos(3)*0.304 -lgd1_pos(4)*0.94 0 lgd1_pos(4)*0.8]; % move the 2nd legend to the desired position
set(lgd3,'box','off','Position',lgd3_pos);
% add text
text(h(1),0.58,2,'open system','FontSize',fs,'Rotation',27)
text(h(1),0.30,13.2,'(unsaturated)','FontSize',fs,'Rotation',53)
text(h(1),0.55,-2.1,'open system','FontSize',fs,'Rotation',20,'color',[0.6 0.6 0.6])
text(h(1),0.26,6.5,'(saturated)','FontSize',fs,'Rotation',45,'color',[0.6 0.6 0.6])
text(h(1),0.75,-10,'closed system','FontSize',fs,'Rotation',7,'Color',[0.3010 0.7450 0.9330])
h(3) = subplot(2,2,3);
% set(h(3),'Position', [0.1 0.1 0.38 0.38]);
hold on
plot(de_ok.f,de_ok.dxsl,'k-','linewidth',2)
plot(de_ok.f,de_ok.dxsv,'k--','linewidth',2)
plot(de_ok.f,de_ok.dxsv_acc,'k:','linewidth',2)
plot(de_oe.f,de_oe.dxsl,'-','color',[0.6 0.6 0.6],'linewidth',2)
plot(de_oe.f,de_oe.dxsv,'--','color',[0.6 0.6 0.6],'linewidth',2)
plot(de_oe.f,de_oe.dxsv_acc,':','color',[0.6 0.6 0.6],'linewidth',2)
plot(de_cl.f,de_cl.dxsl,'-','color',[0.3010 0.7450 0.9330],'linewidth',2)
plot(de_cl.f,de_cl.dxsv,'--','color',[0.3010 0.7450 0.9330],'linewidth',2)
plot(de_cl.f,de_cl.dxsv_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2)
ylim(h(3),[-30 65]);
xlabel('Residual water fraction','fontsize',fs);
ylabel(['\itd / ',char(8240)],'fontsize',fs);
% Rayleigh condensation
h(4) = subplot(2,2,4);
% set(h(4),'Position', [0.56 0.1 0.38 0.38]);
hold on
pc(4) = plot(dc_no.Wconc,dc_no.dxsc,'-','color',[0.6 0.6 0.6],'linewidth',2);
pc(5) = plot(dc_no.Wconc,dc_no.dxs,'--','color',[0.6 0.6 0.6],'linewidth',2);
pc(6) = plot(dc_no.Wconc,dc_no.dxs_acc,':','color',[0.6 0.6 0.6],'linewidth',2);
ylim(gca,[-18 20]);
xlabel('Water mixing ratio / ppmv','fontsize',fs);
set(gca,'XTick',2000:3000:16000);
ylabel(['\itd / ',char(8240)],'fontsize',fs);
% 2nd axis
ax1_pos = get(gca,'position'); % position of first axes
axes('Position',ax1_pos,'XAxisLocation','top','YAxisLocation','right','Color','none','YTickLabel','','FontSize',fs,'TickDir','out');
hold on
pc(1) = plot(dc_ln.f,dc_ln.dxsc,'k-','linewidth',2); % [0.3010 0.7450 0.9330]
pc(2) = plot(dc_ln.f,dc_ln.dxs,'k--','linewidth',2);
pc(3) = plot(dc_ln.f,dc_ln.dxs_acc,'k:','linewidth',2);
% pc(7) = plot(dc_11.f,dc_11.dxsc,'-','color',[0.3010 0.7450 0.9330],'linewidth',2);
% pc(8) = plot(dc_11.f,dc_11.dxs,'--','color',[0.3010 0.7450 0.9330],'linewidth',2);
% pc(9) = plot(dc_11.f,dc_11.dxs_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2);
ylim(gca,[-18 20]);
xlabel('Residual vapour fraction','fontsize',fs);
% 3rd axis
axes('Position',ax1_pos+[0 0 0 -ax1_pos(4)],'XAxisLocation','top','YColor','w','YTick',[],'Color','none','TickDir','out','FontSize',fs);
hold on
pc(7) = plot(dc_ln.Wconc,20*ones(size(dc_ln.Wconc)),'k-'); % [0.3010 0.7450 0.9330]
ylim(gca,[-18 20]);
xlabel('T / °C','fontsize',fs);
set(gca,'XTick',[dc_ln.Wconc(81) dc_ln.Wconc(61:-10:11)],'XTickLabel',[-20 -10:5:15])
h(2) = subplot(2,2,2);
% set(h(2),'Position', [0.56 0.56 0.38 0.38]);
hold on
pc(4) = plot(dc_no.Wconc,dc_no.d18Oc,'-','color',[0.6 0.6 0.6],'linewidth',2); % [0 0.4470 0.7410]
pc(5) = plot(dc_no.Wconc,dc_no.d18O,'--','color',[0.6 0.6 0.6],'linewidth',2);
pc(6) = plot(dc_no.Wconc,dc_no.d18O_acc,':','color',[0.6 0.6 0.6],'linewidth',2);
ylim(h(2),[-60 -2]);
xlabel('Water mixing ratio / ppmv','fontsize',fs);
set(gca,'XTick',2000:3000:16000);
ylabel(['\delta^{18}O / ',char(8240)],'fontsize',fs);
% 2nd axis
ax1_pos = get(gca,'position'); % position of first axes
axes('Position',ax1_pos,'XAxisLocation','top','YAxisLocation','right','Color','none','YTickLabel','','FontSize',fs,'TickDir','out');
hold on
pc(1) = plot(dc_ln.f,dc_ln.d18Oc,'k-','linewidth',2); % [0.3010 0.7450 0.9330]
pc(2) = plot(dc_ln.f,dc_ln.d18O,'k--','linewidth',2);
pc(3) = plot(dc_ln.f,dc_ln.d18O_acc,'k:','linewidth',2);
% pc(7) = plot(dc_11.f,dc_11.d18Oc,'-','color',[0.3010 0.7450 0.9330],'linewidth',2);
% pc(8) = plot(dc_11.f,dc_11.d18O,'--','color',[0.3010 0.7450 0.9330],'linewidth',2);
% pc(9) = plot(dc_11.f,dc_11.d18O_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2);
ylim(gca,[-60 -2]);
xlabel('Residual vapour fraction','fontsize',fs);
lg = legend(h(2),pc(1:3),'instant liquid/ice','residual vapour','total liquid/ice removed');
set(lg,'Position',[0.69 0.66 0.20 0.06],'FontSize',fs);
lg_pos = lg.Position;
% 3rd axis
axes('Position',ax1_pos+[0 -ax1_pos(4)*0.2 0 -ax1_pos(4)*0.8],'XAxisLocation','top','YColor','w','YTick',[],'Color','none');
hold on
pc(7) = plot(dc_ln.Wconc,-2+zeros(size(dc_ln.d18Oc)),'k-'); % [0.3010 0.7450 0.9330]
ylim(gca,[-60 -2]);
xlabel('T / °C','fontsize',fs);
set(gca,'XTick',[dc_ln.Wconc(81) dc_ln.Wconc(61:-10:11)],'XTickLabel',[-20 -10:5:15],'FontSize',fs,'TickDir','out');
lg2 = legend(gca,pc(4:6),'','','','location','east','box','off');
lg2_pos = lg_pos+[-lg_pos(3)*0.405 -lg_pos(4)*0.05 0 lg_pos(4)*0.26]; % move the 2nd legend to the desired position
set(lg2,'box','off','Position',lg2_pos);
% add text
text(h(2),250,-36.9,'super','FontSize',fs,'Rotation',70)
text(h(2),950,-28,'saturation','FontSize',fs,'Rotation',55)
text(h(2),1300,-39,'saturation','FontSize',fs,'Rotation',68,'Color',0.6*[1 1 1])
set([h(1) h(3)],'XDir','reverse','box','on');
set(h(:),'color','none','FontSize',fs,'tickdir','out'); % ,'YAxisLocation','left','ticklength',[0 0],'ygrid','on','xgrid','on'
annotation('textbox',[0.05 0.96 0.02 0.02],'String','(a)','FitBoxToText','on','FontSize',fs+6,'LineStyle','none');
annotation('textbox',[0.49 0.96 0.02 0.02],'String','(b)','FitBoxToText','on','FontSize',fs+6,'LineStyle','none');
annotation('textbox',[0.05 0.48 0.02 0.02],'String','(c)','FitBoxToText','on','FontSize',fs+6,'LineStyle','none');
annotation('textbox',[0.49 0.48 0.02 0.02],'String','(d)','FitBoxToText','on','FontSize',fs+6,'LineStyle','none');
% print('-depsc','-r300','/scratch/home/figs/thesis/Rayleigh_model_plot');
% print('-dpng','-r300','/scratch/home/figs/thesis/Rayleigh_model_plot');