Rayleigh model: Difference between revisions

From gfi
No edit summary
 
(8 intermediate revisions by 2 users not shown)
Line 1: Line 1:




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.  
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 ==
== 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 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 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> 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 formed removal of both phases. First, following the equation above, the total liquid removal is calculated as <math> 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 solid phase. Following the same reasoning as above, the isotope ratio of the accumulated solid removal is calculated by <math> f_iR_\text{i,v} + (1-f_i)R_\text{i,acc} = R_\text{v0} </math>, where <math>R_\text{v0}</math> and <math>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>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> (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 calculated isotopic evolution profiles are shown in Fig. 2.
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> 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> 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> f_iR_\text{i,v} + (1-f_i)R_\text{i,acc} = R_\text{v0} </math>,  
 
where <math>R_\text{v0}</math> and <math>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>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> (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.
 
[[File:Rayleigh_model_box_sketch_v2.png|thumb|800px|'''''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<sub>2</sub><sup>16</sup>O) and ''N<sub>i</sub>'' rare isotopologues (e.g. H<sub>2</sub><sup>18</sup>O or HD<sup>16</sup>O) small amounts of both species, ''dN'' and ''dN<sub>i</sub>'' respectively, are being removed under equilibrium fractionation conditions: <math> 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> 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> R_0N_0 = R(N_0-dN) + αRdN </math>, with the remaining fraction <math> f = (N_0-dN)/N_0 </math>.]]
 
[[File:Rayleigh_model_plot_v2_CG.png|thumb|800px|'''''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 δ<sup>18</sup>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 δ<sup>18</sup>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 δ<sup>18</sup>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<sub>i</sub>''). The supersaturation case (black) takes into account the supersaturation over ice, where the ambient vapour pressure is ''e<sub>v</sub> = S<sub>i</sub>*e<sub>i</sub>''. ''S<sub>i</sub>'' is the defined saturation ratio, as ''e<sub>v</sub>/e<sub>i</sub>'', and here takes the form ''S<sub>i</sub> = 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.]]


{| class="wikitable"
{| class="wikitable"
|+ ''Table 1'' Rayleigh equations to describe isotopic evolution of the remaining reservoir under different conditions.   
|+ ''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>\rho \propto D^{-1} </math>).   
|-
|-
! scope="col"|  
! scope="col"|  
Line 17: Line 35:
! scope="col"| Fractionation factor used
! scope="col"| Fractionation factor used
|-
|-
! scope="row" rowspan="3"| Evaporation
! scope="row" rowspan="4"| Evaporation
| open system saturated
| open system saturated
| <math> R = R_0f^{α-1} </math>
| <math> R = R_0f^{α-1} </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> α = \alpha_l = \frac{R_v}{R_l} <1 </math>
| <math> α = \alpha_{vl}<1 </math>
|-
|-
| open system undersaturated
| open system undersaturated
| <math> R = R_0f^{α-1} </math>
| <math> R = R_0f^{α-1} </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> \alpha = \frac{1}{\alpha_k \alpha_l} </math>, where <math> \alpha_l = \frac{R_l}{R_v} >1 </math>, and <math> \alpha_k = \frac{D}{D'}(1-h)+h > 1 </math>, with <math> h = \frac{e_v}{e_l} </math>   
| <math> \alpha = \frac{1}{\alpha_k \alpha_{lv}} </math>, where <math> \alpha_{lv}>1 </math>, and <math> \alpha_k = \frac{D}{D'}(1-h)+h > 1 </math>, with <math> h = \frac{e_v}{e_l} </math>   
|-
|-
| closed system
| closed system
| <math> R = \frac{R_0}{f+\alpha(1-f)} </math>  
| <math> R = \frac{R_0}{f+\alpha(1-f)} </math>  
| <math> \delta = \frac{\delta_0-(\alpha-1)(1-f)}{f+\alpha(1-f)} </math>  
| <math> \delta = \frac{1+\delta_0}{f+\alpha(1-f)}-1 </math>
| <math> α = \frac{R_v}{R_l} <1 </math>
| <math> α = \alpha_{vl}<1 </math>
|-
| C--G model
| <math> R = R_0f^{\alpha-1} </math>
| <math> \delta = (1+\delta_0)f^{\alpha-1}-1 </math>  
| <math> \alpha = \frac{\alpha_{vl}-h\frac{R_a}{R_l}}{(1-h)\frac{D}{D'}} </math>, where <math> \alpha_{vl} < 1 </math>
|-
|-
! scope="row" rowspan="2"| Condensation
! scope="row" rowspan="2"| Condensation
Line 37: Line 60:
| <math> R = R_0f^{α-1} </math>
| <math> R = R_0f^{α-1} </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> α = \alpha_s = \frac{R_s}{R_v} >1 </math>
| <math> α = \alpha_{sv}>1 </math>
|-
|-
| supersaturation over ice
| supersaturation over ice
| <math> R = R_0f^{α-1} </math>
| <math> R = R_0f^{α-1} </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> δ = (1+δ_0)f^{α-1}-1 </math>
| <math> \alpha = \alpha_k \alpha_s </math>, where <math> \alpha_s = \frac{R_s}{R_v} >1 </math>, and <math> \alpha_k = \frac{S_i}{\alpha_s \frac{D}{D'}(S_i-1)+1} < 1 </math>, with <math> S_i = \frac{e_v}{e_i} </math>   
| <math> \alpha = \alpha_k \alpha_{sv} </math>, where <math> \alpha_{sv}>1 </math>, and <math> \alpha_k = \frac{S_i}{\alpha_s \frac{D}{D'}(S_i-1)+1} < 1 </math>, with <math> S_i = \frac{e_v}{e_i} </math>   
|-
|-
|}
|}


 
----
[[File:Rayleigh_model_box_sketch_v2.png|thumb|800px|'''''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<sub>2</sub><sup>16</sup>O) and ''N<sub>i</sub>'' rare isotopologues (e.g. H<sub>2</sub><sup>18</sup>O or HD<sup>16</sup>O) small amounts of both species, ''dN'' and ''dN<sub>i</sub>'' respectively, are being removed under equilibrium fractionation conditions: <math> 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> 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> R_0N_0 = R(N_0-dN) + αRdN </math>, with the remaining fraction <math> f = (N_0-dN)/N_0 </math>.]]
 
 
[[File:Rayleigh_model_plot_v2.png|thumb|800px|'''''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 δ<sup>18</sup>O = -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 δ<sup>18</sup>O = -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 (''e<sub>i</sub>''). The supersaturation circumstance (black) takes into account the supersaturation over ice where the ambient vapour pressure is ''e<sub>v</sub> = S<sub>i</sub>*e<sub>i</sub>''. ''S<sub>i</sub>'' is the defined saturation ratio, as ''e<sub>v</sub>/e<sub>i</sub>'', and here takes the form ''S<sub>i</sub> = 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 ==
== Model codes in MATLAB ==
Line 65: Line 81:
% simple Rayleigh evaporation model, after Stern and Blisniuk, 2002
% simple Rayleigh evaporation model, after Stern and Blisniuk, 2002
%--------------------------------
%--------------------------------
% input variables: T0 (Kelvin), d18O0 (permil), d2H0 (permil),
% input variables: T0 (Kelvin), d18O0 (permil), d2H0 (permil);
% system (open system = 1, closed system = 0), varargin (relative humidity)
%                 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)
% 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,1)
%              rayleigh_evaporation(20+273.15,-10,-70,0)
%              rayleigh_evaporation(20+273.15,-10,-70,2)
%--------------------------------
%--------------------------------
% HS, Do 12 Dez 2013 15:43:28 CET
% HS, Do 12 Dez 2013 15:43:28 CET
% YW, structured and included kinetic effect, 29 Aug 2020
% 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)
function data = rayleigh_evaporation(T0,d18O0,d2H0,system,h,d18Oa,d2Ha)


% %% dubugging
% %% dubugging
Line 111: Line 129:
R18O0 = (d18O0/1000+1)*R18std;
R18O0 = (d18O0/1000+1)*R18std;
R2H0 = (d2H0/1000+1)*R2std;
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
% % ouput variables
% data.T = [];
% data.T = [];
Line 123: Line 156:
% data.d2Hv_acc = [];
% data.d2Hv_acc = [];
% data.dxsv_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)
%% 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
if system == 1
     fprintf('Rayleigh evaporation process: in an open system (material is immediately removed) \n');
     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('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');
     fprintf('      f    d18Ol    d18Ov    d2Hl    d2Hv    dxsl    dxsv d18Ov_acc dxsv_acc \n');
     n = 1;   % iteration steps
     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
     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
         % modified Rayleigh evaporation model including kinectic effct
         % during the vapour formation at undersaturated conditions (Jouzel and Merlivat 1984)
         % during the vapour formation at undersaturated conditions (Jouzel and Merlivat 1984)
Line 143: Line 218:
         % D_H18O/D_H16O = 0.9723 ± 0.0007
         % D_H18O/D_H16O = 0.9723 ± 0.0007
         % D_HD16O/D_H216O = 0.9755 ± 0.0009
         % D_HD16O/D_H216O = 0.9755 ± 0.0009
         ak_18O = 1/(1/0.9723*(1-h)+h);
        D18 = 0.9723;  % D'/D
         ak_2H = 1/(1/0.9755*(1-h)+h);
        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_18O = a18O(T0) * ak_18O;
         aeff_2H = a2H(T0) * ak_2H;
         aeff_2H = a2H(T0) * ak_2H;
Line 173: Line 250:
%        R18Ov_acc2 = (R18O0-f*R18Ol)/(1-f);
%        R18Ov_acc2 = (R18O0-f*R18Ol)/(1-f);
%        R2Hv_acc2 = (R2H0-f*R2Hl)/(1-f);
%        R2Hv_acc2 = (R2H0-f*R2Hl)/(1-f);
%        data.d18Ov_acc2(n) = (R18Ov_acc2/R18std-1)*1000;
%        data.d18Ov_acc2(i) = (R18Ov_acc2/R18std-1)*1000;
%        data.d2Hv_acc2(n) = (R2Hv_acc2/R2std-1)*1000;
%        data.d2Hv_acc2(i) = (R2Hv_acc2/R2std-1)*1000;
%        data.dxsv_acc2(n) = data.d2Hv_acc2(n)-8*data.d18Ov_acc2(n);
%        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);
         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
         % assign to output
         data.T(n) = T0;
         data.T(i) = T0;
         data.f(n) = f;
         data.f(i) = f;
         data.d18Ol(n) = d18Ol;
        data.R18Ol(i) = R18Ol;
         data.d18Ov(n) = d18Ov;
         data.d18Ol(i) = d18Ol;
         data.d2Hl(n) = d2Hl;
         data.d18Ov(i) = d18Ov;
         data.d2Hv(n) = d2Hv;
        data.R2Hl(i) = R2Hl;
         data.dxsl(n) = dxsl;
         data.d2Hl(i) = d2Hl;
         data.dxsv(n) = dxsv;
         data.d2Hv(i) = d2Hv;
         data.d18Ov_acc(n) = d18Ov_acc;
         data.dxsl(i) = dxsl;
         data.d2Hv_acc(n) = d2Hv_acc;
         data.dxsv(i) = dxsv;
         data.dxsv_acc(n) = dxsv_acc;
         data.d18Ov_acc(i) = d18Ov_acc;
         data.d2Hv_acc(i) = d2Hv_acc;
         data.dxsv_acc(i) = dxsv_acc;


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


%% Rayleigh evaporation process: (2) in closed system (material maintained under two-phase equilibrium)
%% Rayleigh evaporation process: (2) in a "free atmosphere" of fixed humidity and isotope composition (C--G model)
% i.e., the material removed stays in the closed system and stays in equilibrium with the reservoir
if system == 2
% e.g., condensation of vapour to droplets in a cloud (only cloud is formed, but precipitation does not occur)
     fprintf('Rayleigh evaporation process: in an open system (material is immediately removed) \n');
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('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');
     fprintf('      f    d18Ol    d18Ov    d2Hl    d2Hv    dxsl    dxsv d18Ov_acc dxsv_acc \n');
     n = 1;  % iteration steps
    % 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
    % 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
         % isotope ratio of remaining liquid
     n = 1;
         R18Ol = R18O0/(a18O(T0) - f*(a18O(T0)-1)); % Gat1996
    df = 0.02;
        R2Hl = R2H0/(a2H(T0) - f*(a2H(T0)-1));     % Gat1996
    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
         % (1) delta value of remaining liquid
         d18Ol = (R18Ol/R18std-1)*1000;
         d18Ol = (R18Ol/R18std-1)*1000;
Line 215: Line 318:
          
          
         % (2) delta value of the instantaneously removed vapour (in equilibrium with the current remaining liquid)
         % (2) delta value of the instantaneously removed vapour (in equilibrium with the current remaining liquid)
         d18Ov = ((d18Ol/1000+1)*a18O(T0)-1)*1000;
         d18Ov = ((d18Ol/1000+1)*aeff_18O-1)*1000;
         d2Hv = ((d2Hl/1000+1)*a2H(T0)-1)*1000;
         d2Hv = ((d2Hl/1000+1)*aeff_2H-1)*1000;
         dxsv = d2Hv-8*d18Ov;
         dxsv = d2Hv-8*d18Ov;
         % (3) delta value of the accumulated removed vapour
       
         d18Ov_acc = d18Ov;
%         % (3) isotope ratio of the accumulated removed vapour
         d2Hv_acc = d2Hv;
%        R18Ov_acc = R18O0*(1-f^aeff_18O)/(1-f); % Mook 2001, Page 38
         dxsv_acc = dxsv;
%        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);
         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
         % assign to output
         data.T(n) = T0;
         data.T(i) = T0;
         data.f(n) = f;
         data.f(i) = f;
         data.d18Ol(n) = d18Ol;
        data.R18Ol(i) = R18Ol;
         data.d18Ov(n) = d18Ov;
         data.d18Ol(i) = d18Ol;
         data.d2Hl(n) = d2Hl;
         data.d18Ov(i) = d18Ov;
         data.d2Hv(n) = d2Hv;
        data.R2Hl(i) = R2Hl;
         data.dxsl(n) = dxsl;
         data.d2Hl(i) = d2Hl;
         data.dxsv(n) = dxsv;
         data.d2Hv(i) = d2Hv;
         data.d18Ov_acc(n) = d18Ov_acc;
         data.dxsl(i) = dxsl;
         data.d2Hv_acc(n) = d2Hv_acc;
         data.dxsv(i) = dxsv;
         data.dxsv_acc(n) = dxsv_acc;
         data.d18Ov_acc(i) = d18Ov_acc;
         data.d2Hv_acc(i) = d2Hv_acc;
         data.dxsv_acc(i) = dxsv_acc;
          
          
         % iterate to next step
         % iterate to next step
         n = n+1;
         i = i+1;
       
     end
     end
end
end
Line 255: Line 374:
% Description:
% Description:
% - simple Rayleigh condensation model, after Stern and Blisniuk, 2002
% - simple Rayleigh condensation model, after Stern and Blisniuk, 2002
% - temperature dependent fractionation factor varies during cooling process
% - temperature dependent fractionation factor varies at each step during cooling process
%  (1) isobaric cooling: e.g., air parcel advects to cold regions
%  (1) isobaric cooling: e.g., air parcel advects to cold regions
%  (2) or pseudoadiabatic cooling: e.g., updraft or lift of surface air parcel
%  (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),
% input variables: T0 (Kelvin), RH0 (%), d18O0 (permil), d2H0 (permil),
Line 413: Line 533:
     data.rh(n) = rh;
     data.rh(n) = rh;
     data.Wconc(n) = wconc;
     data.Wconc(n) = wconc;
    data.R18O(n) = R18O0;
    data.R2H(n) = R2H0;
     data.d18O(n) = d18O0;
     data.d18O(n) = d18O0;
     data.dD(n) = d2H0;
     data.dD(n) = d2H0;
Line 438: Line 560:
     f = ws/ws0;
     f = ws/ws0;
      
      
     % isotope ratio of remaining vapour
     if n == 1
    R18Ov = R18O0*f^(a18O(T)-1);
        f1 = ws/ws0;
     R2Hv = R2H0*f^(a2H(T)-1);
        % 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
     % (1) delta value of remaining vapour
     d18Ov = (R18Ov/R18std-1)*1000;
     d18Ov = (R18Ov/R18std-1)*1000;
Line 467: Line 598:
     data.rh(n) = rh;
     data.rh(n) = rh;
     data.Wconc(n) = wconc;
     data.Wconc(n) = wconc;
    data.R18O(n) = R18Ov;
    data.R2H(n) = R2Hv;
     data.d18O(n) = d18Ov;
     data.d18O(n) = d18Ov;
     data.dD(n) = d2Hv;
     data.dD(n) = d2Hv;
Line 504: Line 637:
     f = ws/ws0;      % the total fraction of remaining vapour through all stages
     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
     fi = ws/ws0_ice; % the fraction of remaing vapour in the ice condensation stage only
   
       
     % modified Rayleigh condensation model including kinectic effct
     % modified Rayleigh condensation model including kinectic effct
     % during the ice formation at supersaturated conditions (Jouzel and Merlivat 1984)
     % during the ice formation at supersaturated conditions (Jouzel and Merlivat 1984)
Line 515: Line 648:
     aeff_2H = a2H(T) * ak_2H;
     aeff_2H = a2H(T) * ak_2H;
      
      
     % isotope ratio of remaining vapour
     if n == 1
    R18Ov = R18O0*fi^(aeff_18O-1); % Jouzel and Merlivat, 1984
        f1 = ws/ws0;
    R2Hv = R2H0*fi^(aeff_2H-1);   % Jouzel and Merlivat, 1984
        % 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
     % (1) delta value of remaining vapour
     d18Ov = (R18Ov/R18std-1)*1000;
     d18Ov = (R18Ov/R18std-1)*1000;
Line 531: Line 673:
     % (i) mass conservation of rare isotopes between the reservoir and formed removal:
     % (i) mass conservation of rare isotopes between the reservoir and formed removal:
     % fi*Rv + (1-fi)*Ri_acc = Rv0,
     % fi*Rv + (1-fi)*Ri_acc = Rv0,
     % where fi is the fraction of remaing vapour in the ice condensation stage, and
     % 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
     % Rv0 is the isotope ratio of the remaining vapour after liquid condensation stage
     % after rearranging, we have: Ri_acc = (Rv0-fi*Rv)/(1-fi)
     % after rearranging, we have: Ri_acc = (Rv0-fi*Rv)/(1-fi)
Line 554: Line 696:
     data.rh(n) = rh;
     data.rh(n) = rh;
     data.Wconc(n) = wconc;
     data.Wconc(n) = wconc;
    data.R18O(n) = R18Ov;
    data.R2H(n) = R2Hv;
     data.d18O(n) = d18Ov;
     data.d18O(n) = d18Ov;
     data.dD(n) = d2Hv;
     data.dD(n) = d2Hv;
Line 609: Line 753:
de_oe = rayleigh_evaporation(T0,d18O0,d2H0,1);  % in an open system, saturated (equilibrium)
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_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
% Rayleigh condensation
T0 = 20+273.15;
T0 = 20+273.15;
Line 651: Line 797:
pe(8) = plot(de_cl.f,de_cl.d18Ov,'--','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(9) = plot(de_cl.f,de_cl.d18Ov_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2);
ylim(h(1),[-27 40]);
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);
xlabel('Residual water fraction','fontsize',fs);
ylabel(['\delta^{18}O / ',char(8240)],'fontsize',fs);
ylabel(['\delta^{18}O / ',char(8240)],'fontsize',fs);
Line 657: Line 806:
lgd1_pos = lgd1.Position;
lgd1_pos = lgd1.Position;
lgd1.FontSize = fs;
lgd1.FontSize = fs;
% lgd1.Visible = 'off';
% 2nd axes
% 2nd axes
a2 = axes('position',get(h(1),'position'),'visible','off'); % create the same axes but set it invisible
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 = 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
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);
set(lgd2,'box','off','Position',lgd2_pos);
% 3rd axes
% 3rd axes
a3 = axes('position',get(gca,'position'),'visible','off'); % create the same axes but set it invisible
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 = 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
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);
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
% add text
text(h(1),0.58,2,'open system','FontSize',fs,'Rotation',27)
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.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.60,-3,'open','FontSize',fs,'Rotation',15,'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.48,-0.5,'(saturated)','FontSize',fs,'Rotation',20,'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])
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);
h(3) = subplot(2,2,3);
Line 686: Line 844:
plot(de_cl.f,de_cl.dxsv,'--','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_cl.f,de_cl.dxsv_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2)
ylim(h(3),[-30 65]);
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);
xlabel('Residual water fraction','fontsize',fs);
ylabel(['\itd / ',char(8240)],'fontsize',fs);
ylabel(['\itd / ',char(8240)],'fontsize',fs);
Line 697: Line 859:
pc(5) = plot(dc_no.Wconc,dc_no.dxs,'--','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);   
pc(6) = plot(dc_no.Wconc,dc_no.dxs_acc,':','color',[0.6 0.6 0.6],'linewidth',2);   
ylim(gca,[-18 20]);
ylim(gca,[-20 50]);
xlabel('Water mixing ratio / ppmv','fontsize',fs);
xlabel('Water mixing ratio / ppmv','fontsize',fs);
set(gca,'XTick',2000:3000:16000);
set(gca,'XTick',2000:3000:16000);
Line 711: Line 873:
% pc(8) = plot(dc_11.f,dc_11.dxs,'--','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);
% pc(9) = plot(dc_11.f,dc_11.dxs_acc,':','color',[0.3010 0.7450 0.9330],'linewidth',2);
ylim(gca,[-18 20]);
ylim(gca,[-20 50]);
xlabel('Residual vapour fraction','fontsize',fs);
xlabel('Residual vapour fraction','fontsize',fs);
% 3rd axis
% 3rd axis
Line 717: Line 879:
hold on
hold on
pc(7) = plot(dc_ln.Wconc,20*ones(size(dc_ln.Wconc)),'k-'); % [0.3010 0.7450 0.9330]
pc(7) = plot(dc_ln.Wconc,20*ones(size(dc_ln.Wconc)),'k-'); % [0.3010 0.7450 0.9330]
ylim(gca,[-18 20]);
ylim(gca,[-20 50]);
xlabel('T / °C','fontsize',fs);
xlabel('T / °C','fontsize',fs);
set(gca,'XTick',[dc_ln.Wconc(81) dc_ln.Wconc(61:-10:11)],'XTickLabel',[-20 -10:5:15])
set(gca,'XTick',[dc_ln.Wconc(81) dc_ln.Wconc(61:-10:11)],'XTickLabel',[-20 -10:5:15])
Line 744: Line 906:
xlabel('Residual vapour fraction','fontsize',fs);
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(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);
set(lg,'Position',[0.69 0.66 0.20 0.06],'FontSize',fs);
lg_pos = lg.Position;
lg_pos = lg.Position;
Line 754: Line 917:
set(gca,'XTick',[dc_ln.Wconc(81) dc_ln.Wconc(61:-10:11)],'XTickLabel',[-20 -10:5:15],'FontSize',fs,'TickDir','out');
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 = 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
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);
set(lg2,'box','off','Position',lg2_pos);
% add text
% add text
text(h(2),250,-36.9,'super','FontSize',fs,'Rotation',70)
text(h(2),1300,-36,'super','FontSize',fs,'Rotation',65)
text(h(2),950,-28,'saturation','FontSize',fs,'Rotation',55)
text(h(2),2150,-27.5,'saturation','FontSize',fs,'Rotation',50)
text(h(2),1300,-39,'saturation','FontSize',fs,'Rotation',68,'Color',0.6*[1 1 1])
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(1) h(3)],'XDir','reverse','box','on');

Latest revision as of 20:26, 26 January 2021


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. 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 here 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). 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 δ18O = -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 δ18O = -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 (ei). The supersaturation case (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 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');
Retrieved from "http://gfi.wiki.uib.no/index.php?title=Rayleigh_model&oldid=1448"