Rayleigh model

The 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 fractionation process.

Model description

The system applying the Rayleigh equation is often an open system, from which the newly formed phase is immediately removed. For example, the evaporation of natural water bodies, and the formation of precipitation that immediately leaves the system, can be regarded as examples for open systems. However, with some modifications, Rayleigh models 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 for this would be the condensation of vapour to droplets in a cloud, without formation and removal of precipitation.

The isotopic enrichment or depletion by the Rayleigh process for both open and closed systems can be mathematically obtained by different approaches (Fig. 1). The Rayleigh equations that are used to describe the 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 remaining in the reservoir, and the removed component:

[math]\displaystyle{ fR_\text{remain} + (1-f)R_\text{remove,acc} = R_0 }[/math].

In the case of condensation involving both liquid and solid phases, the accumulation is a combination of the removed mass fraction of both phases. First, following the equation above, the total liquid removal is calculated as

[math]\displaystyle{ R_\text{l,acc} = (R_0-fR_v)/(1-f) }[/math].

Then, we consider a separated Rayleigh process by only taking into account of the ongoing formation of the 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 the 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]. Corresponding examples of the calculated isotopic evolution profiles are shown in Fig. 2.

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. H₂¹⁶O) and N_i rare isotopologues (e.g. H₂¹⁸O or HD¹⁶O) small amounts of both species, dN and dN_i 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 here at 20 °C (thus a constant fractionation factor) for initial liquid compositions of δ¹⁸O = -10 ‰ and δD = -70 ‰ (d = 10 ‰), in open system with unsaturated (RH = 75 %; black) or saturated (grey) environment, and closed system (blue). In addition, the result of C--G evaporation model (i.e. open system including feedback from a "free atmosphere" which has fixed isotope compositions of δ¹⁸O = -13 ‰ and δD = -94 ‰ (d = 10 ‰) and a fixed relative humidity (RH = 75 %)) is also presented (red). Rayleigh condensation occurs here under continuous cooling (thus also a contentiously changing fractionation factor) from T = 20 °C and RH = 75 %. The initial vapour compositions are δ¹⁸O = -13 ‰ and δD = -94 ‰ (d = 10 ‰). For the condensation to ice below 0 °C, two cases are presented. The saturation case (grey) is a classical Rayleigh process where ice crystals from from vapour under equilibrium conditions, using the saturation pressure over ice (e_i). The supersaturation case (black) takes into account the supersaturation over ice, where the ambient vapour pressure is e_v = S_i*e_i. S_i is the defined saturation ratio, as e_v/e_i, and here takes the form S_i = 1-0.004T after Risi et al (2010a). In this case, the fractionation factor combining equilibrium and kinetic effects given by Jouzel & Merlivat (1984) is used. The equations used in this figure can be found in Table 1.

Table 1 Rayleigh equations to calculate isotopic evolution of the remaining reservoir under different conditions. The transport resistance is here assumed to be purely due to diffusion (i.e. [math]\displaystyle{ \rho \propto D^{-1} }[/math]).

	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_{vl}\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_{lv}} }[/math], where [math]\displaystyle{ \alpha_{lv}\gt 1 }[/math], and [math]\displaystyle{ \alpha_k = \frac{D}{D'}(1-h)+h \gt 1 }[/math], with [math]\displaystyle{ h = \frac{e_v}{e_l} }[/math]
	closed system	[math]\displaystyle{ R = \frac{R_0}{f+\alpha(1-f)} }[/math]	[math]\displaystyle{ \delta = \frac{1+\delta_0}{f+\alpha(1-f)}-1 }[/math]	[math]\displaystyle{ α = \alpha_{vl}\lt 1 }[/math]
	C--G model	[math]\displaystyle{ R = R_0f^{\alpha-1} }[/math]	[math]\displaystyle{ \delta = (1+\delta_0)f^{\alpha-1}-1 }[/math]	[math]\displaystyle{ \alpha = \frac{\alpha_{vl}-h\frac{R_a}{R_l}}{(1-h)\frac{D}{D'}} }[/math], where [math]\displaystyle{ \alpha_{vl} \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_{sv}\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_{sv} }[/math], where [math]\displaystyle{ \alpha_{sv}\gt 1 }[/math], and [math]\displaystyle{ \alpha_k = \frac{S_i}{\alpha_s \frac{D}{D'}(S_i-1)+1} \lt 1 }[/math], with [math]\displaystyle{ S_i = \frac{e_v}{e_i} }[/math]

---

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: closed system = 0, open system = 1, C--G model = 2; 
%                  optional: relative humidity (default = 1), d18Oa (permil), d2Ha (permil)
% 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,2)
%--------------------------------
% HS, Do 12 Dez 2013 15:43:28 CET
% YW, structured and included the case of kinetic effect 
% and the case of C--G model, 29 Aug 2020
%--------------------------------

function data = rayleigh_evaporation(T0,d18O0,d2H0,system,h,d18Oa,d2Ha)

% %% 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;
% relative humidity
if ~exist('h','var')
    h = 1;
end
% isotope composition of "free atmosphere"
% if not existing, default it to the equilibrium vapour of global average precipitation 
if ~exist('d18Oa','var')
    d18Oa = -13;  
end
if ~exist('d2Ha','var')
    d2Ha = -94;  
end
R18Oa = (d18Oa/1000+1)*R18std;
R2Ha = (d2Ha/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: (0) 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
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');
    i = 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(i) = T0;
        data.f(i) = f;
        data.R18Ol(i) = R18Ol;
        data.d18Ol(i) = d18Ol;
        data.d18Ov(i) = d18Ov;
        data.R2Hl(i) = R2Hl;
        data.d2Hl(i) = d2Hl;
        data.d2Hv(i) = d2Hv;
        data.dxsl(i) = dxsl;
        data.dxsv(i) = dxsv;
        data.d18Ov_acc(i) = d18Ov_acc;
        data.d2Hv_acc(i) = d2Hv_acc;
        data.dxsv_acc(i) = dxsv_acc;
        
        % iterate to next step
        i = i+1;
    end
end

%% Rayleigh evaporation process: (1) in an open system (material is immediately removed)
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');
    i = 1; % iteration steps
    n = 1; % turbulent exponent in D^n, [0 1], controls the ratio of diffusive transport over turbulent transport of water molecules (Pfahl and Wernli, 2009)
    for f = 1:-0.02:0 % remaining liquid fraction
       
        % 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
        D18 = 0.9723;  % D'/D
        D2 = 0.9755;   % D'/D
        ak_18O = 1/(1/D18^n*(1-h)+h);
        ak_2H = 1/(1/D2^n*(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(i) = (R18Ov_acc2/R18std-1)*1000;
%         data.d2Hv_acc2(i) = (R2Hv_acc2/R2std-1)*1000;
%         data.dxsv_acc2(i) = data.d2Hv_acc2(i)-8*data.d18Ov_acc2(i);
        
        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(i) = T0;
        data.f(i) = f;
        data.R18Ol(i) = R18Ol;
        data.d18Ol(i) = d18Ol;
        data.d18Ov(i) = d18Ov;
        data.R2Hl(i) = R2Hl;
        data.d2Hl(i) = d2Hl;
        data.d2Hv(i) = d2Hv;
        data.dxsl(i) = dxsl;
        data.dxsv(i) = dxsv;
        data.d18Ov_acc(i) = d18Ov_acc;
        data.d2Hv_acc(i) = d2Hv_acc;
        data.dxsv_acc(i) = dxsv_acc;

        % iterate to next step
        i = i+1;
    end
end

%% Rayleigh evaporation process: (2) in a "free atmosphere" of fixed humidity and isotope composition (C--G model)
if system == 2
    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');
    % turbulent exponent in D^n, [0 1], controls the ratio of diffusive transport over turbulent transport of water molecules (Pfahl and Wernli, 2009)
    % n~0.5 for an open water body under natural conditions, n~0.58 for a falling raindrop (Stewart1975), n~1 for soils and leaves 
    n = 1; 
    df = 0.02;
    i = 1;   % iteration steps
    
    for f = 1:-df:0 % remaining liquid fraction

        % 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
        D18 = 0.9723;  % D'/D
        D2 = 0.9755;   % D'/D
        if i==1
            aeff_18O = (a18O(T0) - h*R18Oa/R18O0) / ((1-h)*1/D18^n);
            aeff_2H = (a2H(T0) - h*R2Ha/R2H0) / ((1-h)*1/D2^n);
            % isotope ratio of remaining liquid
            f1 = 1;
            R18Ol = R18O0*f1^(aeff_18O-1);
            R2Hl = R2H0*f1^(aeff_2H-1);
        else
            aeff_18O = (a18O(T0) - h*R18Oa/data.R18Ol(i-1)) / ((1-h)*1/D18^n);
            aeff_2H = (a2H(T0) - h*R2Ha/data.R2Hl(i-1)) / ((1-h)*1/D2^n);
            % isotope ratio of remaining liquid
            f1 = 1-df; % fraction remaining relative to previous stage
            R18Ol = data.R18Ol(i-1)*f1^(aeff_18O-1);
            R2Hl = data.R2Hl(i-1)*f1^(aeff_2H-1);
        end
       
        % (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_acc = (R18O0-f*R18Ol)/(1-f);
        R2Hv_acc = (R2H0-f*R2Hl)/(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;
        
        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(i) = T0;
        data.f(i) = f;
        data.R18Ol(i) = R18Ol;
        data.d18Ol(i) = d18Ol;
        data.d18Ov(i) = d18Ov;
        data.R2Hl(i) = R2Hl;
        data.d2Hl(i) = d2Hl;
        data.d2Hv(i) = d2Hv;
        data.dxsl(i) = dxsl;
        data.dxsv(i) = dxsv;
        data.d18Ov_acc(i) = d18Ov_acc;
        data.d2Hv_acc(i) = d2Hv_acc;
        data.dxsv_acc(i) = dxsv_acc;
        
        % iterate to next step
        i = i+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 at each step 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
% - therefore, the isotope ratio is calculated as: R(n+1) = R(n)(q(n+1)/q(n))^(a-1)
%--------------------------------
% 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.R18O(n) = R18O0;
    data.R2H(n) = R2H0;
    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;
    
    if n == 1
        f1 = ws/ws0;
        % isotope ratio of remaining vapour
        R18Ov = R18O0*f1^(a18O(T)-1); % R(n) = R(n-1)(q(n)/q(n-1))^(a-1)
        R2Hv = R2H0*f1^(a2H(T)-1);
    else
        f1 = wconc/data.Wconc(n-1); % q(n)/q(n-1)
        % isotope ratio of remaining vapour
        R18Ov = data.R18O(n-1)*f1^(a18O(T)-1);  % R(n) = R(n-1)(q(n)/q(n-1))^(a-1)
        R2Hv = data.R2H(n-1)*f1^(a2H(T)-1);
    end
    
    % (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.R18O(n) = R18Ov;
    data.R2H(n) = R2Hv;
    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;
    
    if n == 1
        f1 = ws/ws0;
        % isotope ratio of remaining vapour
        R18Ov = R18O0*f1^(aeff_18O-1); % R(n) = R(n-1)(q(n)/q(n-1))^(a-1)
        R2Hv = R2H0*f1^(aeff_2H-1);
    else
        f1 = wconc/data.Wconc(n-1); % q(n)/q(n-1)
        % isotope ratio of remaining vapour
        R18Ov = data.R18O(n-1)*f1^(aeff_18O-1); % R(n) = R(n-1)(q(n)/q(n-1))^(a-1), and Jouzel and Merlivat, 1984
        R2Hv = data.R2H(n-1)*f1^(aeff_2H-1);
    end
    
    % (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 remaining 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.R18O(n) = R18Ov;
    data.R2H(n) = R2Hv;
    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
de_CG = rayleigh_evaporation(T0,d18O0,d2H0,2,h); % C--G model

% 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);
pe(10) = plot(de_CG.f,de_CG.d18Ol,'-','color',[0.8500 0.3250 0.0980],'linewidth',2); 
pe(11) = plot(de_CG.f,de_CG.d18Ov,'--','color',[0.8500 0.3250 0.0980],'linewidth',2);
pe(12) = plot(de_CG.f,de_CG.d18Ov_acc,':','color',[0.8500 0.3250 0.0980],'linewidth',2);
ylim(h(1),[-30 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;
% lgd1.Visible = 'off';
% 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),'residual water','instant vapour','total vapour removed','Location','northwest');
% lgd2.FontSize = fs;
lgd2 = legend(a2,pe(4:6),'','','','location','northwest','box','off');
lgd2_pos = lgd1_pos+[-lgd1_pos(3)*0.304 -lgd1_pos(4)*0.58 0 lgd1_pos(4)*0.5]; % 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.79 0 lgd1_pos(4)*0.5]; % move the 2nd legend to the desired position
set(lgd3,'box','off','Position',lgd3_pos);
% % 4th axes
a4 = axes('position',get(gca,'position'),'visible','off'); % create the same axes but set it invisible
lgd4 = legend(a4,pe(10:12),'','','','location','northwest','box','off');
lgd4_pos = lgd1_pos+[-lgd1_pos(3)*0.304 -lgd1_pos(4)*0.88 0 lgd1_pos(4)*0.5]; % move the 2nd legend to the desired position
set(lgd4,'box','off','Position',lgd4_pos);
% add text
text(h(1),0.41,8.5,'open','FontSize',fs,'Rotation',30)
text(h(1),0.30,13.2,'(unsaturated)','FontSize',fs,'Rotation',53)
text(h(1),0.60,-3,'open','FontSize',fs,'Rotation',15,'color',[0.6 0.6 0.6])
text(h(1),0.48,-0.5,'(saturated)','FontSize',fs,'Rotation',20,'color',[0.6 0.6 0.6])
text(h(1),0.75,-10,'closed','FontSize',fs,'Rotation',7,'Color',[0.3010 0.7450 0.9330])
text(h(1),0.8,1.5,'C-G model','color',[0.8500 0.3250 0.0980],'FontSize',fs,'Rotation',15)

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)
plot(de_CG.f,de_CG.dxsl,'-','color',[0.8500 0.3250 0.0980],'linewidth',2)
plot(de_CG.f,de_CG.dxsv,'--','color',[0.8500 0.3250 0.0980],'linewidth',2)
plot(de_CG.f,de_CG.dxsv_acc,':','color',[0.8500 0.3250 0.0980],'linewidth',2)
ylim(h(3),[-60 70]);
set(h(3),'ytick',-60:20: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,[-20 50]);
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,[-20 50]);
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,[-20 50]);
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');
% lg = legend(h(2),pc(4:6),'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.02 0 lg_pos(4)*0.117]; % move the 2nd legend to the desired position
set(lg2,'box','off','Position',lg2_pos);
% add text
text(h(2),1300,-36,'super','FontSize',fs,'Rotation',65)
text(h(2),2150,-27.5,'saturation','FontSize',fs,'Rotation',50)
text(h(2),1200,-25,'saturation','FontSize',fs,'Rotation',50,'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');