Phase Change

function [dummy] = neumann

% --------------------------------------------------------
% Analytical solution for Neumann problem (Lunardini 1981)
% --------------------------------------------------------

% soil depths (meters)

nsoi = 4;
depth(1) = 0.25;
depth(2) = 0.55;
depth(3) = 0.85;
depth(4) = 1.15;

% number of days to simulate

ndays = 60;

% time-invariant surface temperature (deg C)

ts = -10;

% initial soil temperature (deg C)

for i = 1:nsoi
   t0(i) = 2;
end

% freezing temperature (deg C)

tf = 0;

% soil dependent constants
% a1 - frozen diffusivity (cm**2/hr -> m**2/s)
% a2 - unfrozen diffusivity (cm**2/hr -> m**2/s)
% rm - "m" from Jumikis (1966) and adjust for above unit conversion
% rg - "gamma" from Lunardini (1981)

a1 = 42.55 / 3600 / (100*100);
a2 = 23.39 / 3600 / (100*100);
rm = 3.6 / (60 * 100);
rg = rm / (2 * sqrt(a1));

% save output for day zero (i.e., initial conditions)

m = 1;
iday_out(m) = 0;
xfr_out(m) = 0;
t1_out(m) = t0(1);
t2_out(m) = t0(2);
t3_out(m) = t0(3);
t4_out(m) = t0(4);
z1_out(m) = depth(1) * 100;
z2_out(m) = depth(2) * 100;
z3_out(m) = depth(3) * 100;
z4_out(m) = depth(4) * 100;

% begin time loop

for iday = 1:ndays

   % time (seconds)

   time = iday * 24 * 3600;

   % depth of frost penetration XF (meters) based on time (seconds)

   xf = 2 * rg * sqrt(a1*time);

   % calculate soil temperatures given XF and time

   for i = 1:nsoi
      if (depth(i) <= xf)
         x = depth(i) / (2*sqrt(a1*time));
         y = rg;
         t(i) = ts + (tf - ts) * erf(x) / erf(y);
      else
         x = depth(i) / (2*sqrt(a2*time));
         y = rg * sqrt(a1/a2);
         t(i) = t0(i) - (t0(i) - tf) * (1-erf(x)) / (1-erf(y));
      end
   end

   % save for output

   m = m + 1;
   iday_out(m) = iday;
   xf_out(m) = xf * 100;
   t1_out(m) = t(1);
   t2_out(m) = t(2);
   t3_out(m) = t(3);
   t4_out(m) = t(4);
   z1_out(m) = depth(1) * 100;
   z2_out(m) = depth(2) * 100;
   z3_out(m) = depth(3) * 100;
   z4_out(m) = depth(4) * 100;

end

A = [iday_out; xf_out; z1_out; t1_out; z2_out; t2_out; z3_out; t3_out; z4_out; t4_out];
fileID = fopen('data_analytical.txt','w');
fprintf(fileID,'%10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n','day','z0C','z1','t1','z2','t2','z3','t3','z4','t4');
fprintf(fileID,'%10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n', A);
fclose(fileID);

dummy = 0;

function [soilvar] = phase_change (physcon, soilvar, dt)

% Adjust temperatures for phase change. Freeze or melt ice using
% energy excess or deficit needed to change temperature to the
% freezing point.

% ------------------------------------------------------
% Input
%   dt                      ! Time step (s)
%   physcon.hfus            ! Heat of fusion for water at 0 C (J/kg)
%   physcon.tfrz            ! Freezing point of water (K)
%   soilvar.nsoi            ! Number of soil layers
%   soilvar.dz              ! Soil layer thickness (m)
%   soilvar.cv              ! Volumetric heat capacity (J/m3/K)
%
% Input/output
%   soilvar.tsoi            ! Soil temperature (K)
%   soilvar.h2osoi_liq      ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice      ! Frozen water, ice (kg H2O/m2)
%
% Output
%   soilvar.hfsoi           ! Soil phase change energy flux (W/m2)
% ------------------------------------------------------

% --- Initialize total soil heat of fusion to zero

soilvar.hfsoi = 0;

% --- Now loop over all soil layers to calculate phase change

for i = 1:soilvar.nsoi

   % --- Save variables prior to phase change

   wliq0 = soilvar.h2osoi_liq(i);     % Amount of liquid water before phase change
   wice0 = soilvar.h2osoi_ice(i);     % Amount of ice before phase change
   wmass0 = wliq0 + wice0;            % Amount of total water before phase change
   tsoi0 = soilvar.tsoi(i);           % Soil temperature before phase change

   % --- Identify melting or freezing layers and set temperature to freezing

   % Default condition is no phase change (imelt = 0)

   imelt = 0;

   % Melting: if ice exists above melt point, melt some to liquid.
   % Identify melting by imelt = 1

   if (soilvar.h2osoi_ice(i) > 0 & soilvar.tsoi(i) > physcon.tfrz)
      imelt = 1;
      soilvar.tsoi(i) = physcon.tfrz;
   end

   % Freezing: if liquid exists below melt point, freeze some to ice.
   % Identify freezing by imelt = 2

   if (soilvar.h2osoi_liq(i) > 0 & soilvar.tsoi(i) < physcon.tfrz)
      imelt = 2;
      soilvar.tsoi(i) = physcon.tfrz;
   end

   % --- Calculate energy for freezing or melting

   % The energy for freezing or melting (W/m2) is assessed from the energy
   % excess or deficit needed to change temperature to the freezing point.
   % This is a potential energy flux, because cannot melt more ice than is
   % present or freeze more liquid water than is present.
   %
   % heat_flux_pot > 0: freezing; heat_flux_pot < 0: melting

   if (imelt > 0)
      heat_flux_pot = (soilvar.tsoi(i) - tsoi0) * soilvar.cv(i) * soilvar.dz(i) / dt;
   else
      heat_flux_pot = 0;
   end

   % Maximum energy for melting or freezing (W/m2)

   if (imelt == 1)
      heat_flux_max = -soilvar.h2osoi_ice(i) * physcon.hfus / dt;
   end

   if (imelt == 2)
      heat_flux_max = soilvar.h2osoi_liq(i) * physcon.hfus / dt;
   end

   % --- Now freeze or melt ice

   if (imelt > 0)

      % Change in ice (kg H2O/m2/s): freeze (+) or melt (-)

      ice_flux = heat_flux_pot / physcon.hfus;

      % Update ice (kg H2O/m2)

      soilvar.h2osoi_ice(i) = wice0 + ice_flux * dt;

      % Cannot melt more ice than is present

      soilvar.h2osoi_ice(i) = max(0, soilvar.h2osoi_ice(i));

      % Ice cannot exceed total water that is present

      soilvar.h2osoi_ice(i) = min(wmass0, soilvar.h2osoi_ice(i));

      % Update liquid water (kg H2O/m2) for change in ice

      soilvar.h2osoi_liq(i) = max(0, (wmass0-soilvar.h2osoi_ice(i)));

      % Actual energy flux from phase change (W/m2). This is equal to
      % heat_flux_pot except if tried to melt too much ice.

      heat_flux = physcon.hfus * (soilvar.h2osoi_ice(i) - wice0) / dt;

      % Sum energy flux from phase change (W/m2)

      soilvar.hfsoi = soilvar.hfsoi + heat_flux;

      % Residual energy not used in phase change is added to soil temperature

      residual = heat_flux_pot - heat_flux;
      soilvar.tsoi(i) = soilvar.tsoi(i) - residual * dt / (soilvar.cv(i) * soilvar.dz(i));

      % Error check: make sure actual phase change does not exceed permissible phase change

      if (abs(heat_flux) > abs(heat_flux_max))
         error ('Soil temperature energy conservation error: phase change')
      end

      % Freezing: make sure actual phase change does not exceed permissible phase change
      % and that the change in ice does not exceed permissible change

      if (imelt == 2)

         % Energy flux (W/m2)

         constraint = min(heat_flux_pot, heat_flux_max);
         err = heat_flux - constraint;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: freezing energy flux')
         end

         % Change in ice (kg H2O/m2)

         err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: freezing ice flux')
         end
      end

      % Thawing: make sure actual phase change does not exceed permissible phase change
      % and that the change in ice does not exceed permissible change

      if (imelt == 1)

         % Energy flux (W/m2)

         constraint = max(heat_flux_pot, heat_flux_max);
         err = heat_flux - constraint;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: thawing energy flux')
         end

         % Change in ice (kg H2O/m2)

         err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
         if (abs(err) > 1e-03)
            error ('Soil temperature energy conservation error: thawing ice flux')
         end
      end

   end

end

function [soilvar] = soil_temperature (physcon, soilvar, tsurf, dt)

% Use an implicit formulation with the surface boundary condition specified
% as the surface temperature to solve for soil temperatures at time n+1.
%
% Calculate soil temperatures as:
%
%      dT   d     dT 
%   cv -- = -- (k --)
%      dt   dz    dz 
%
% where: T = temperature (K)
%        t = time (s)
%        z = depth (m)
%        cv = volumetric heat capacity (J/m3/K)
%        k = thermal conductivity (W/m/K)
%
% Set up a tridiagonal system of equations to solve for T at time n+1, 
% where the temperature equation for layer i is
%
%   d_i = a_i [T_i-1] n+1 + b_i [T_i] n+1 + c_i [T_i+1] n+1
%
% For soil layers undergoing phase change, set T_i = Tf (freezing) and use
% excess energy to freeze or melt ice:
%
%   Hf_i = (Tf - [T_i] n+1) * cv_i * dz_i / dt
%
% During the phase change, the unfrozen and frozen soil water
% (h2osoi_liq, h2osoi_ice) are adjusted.
%
% Or alternatively, use the apparent heat capacity method to
% account for phase change. In this approach, h2osoi_liq
% and h2osoi_ice are not calculated.
%
% ------------------------------------------------------
% Input
%   tsurf                   ! Surface temperature (K)
%   dt                      ! Time step (s)
%   soilvar.method          ! Use excess heat or apparent heat capacity for phase change
%   soilvar.nsoi            ! Number of soil layers
%   soilvar.z               ! Soil depth (m)
%   soilvar.z_plus_onehalf  ! Soil depth (m) at i+1/2 interface between layers i and i+1
%   soilvar.dz              ! Soil layer thickness (m)
%   soilvar.dz_plus_onehalf ! Thickness (m) between between i and i+1
%   soilvar.tk              ! Thermal conductivity (W/m/K)
%   soilvar.cv              ! Heat capacity (J/m3/K)
%
% Input/output
%   soilvar.tsoi            ! Soil temperature (K)
%   soilvar.h2osoi_liq      ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice      ! Frozen water, ice (kg H2O/m2)
%
% Output
%   soilvar.gsoi            ! Energy flux into soil (W/m2)
%   soilvar.hfsoi           ! Soil phase change energy flux (W/m2)
% ------------------------------------------------------

% solution = 'Crank-Nicolson'; % Use Crank-Nicolson solution
solution = 'implicit';       % Use implicit solution

% --- Save current soil temperature for energy conservation check

for i = 1:soilvar.nsoi
   tsoi0(i) = soilvar.tsoi(i);
end

% --- Thermal conductivity at interface (W/m/K)

for i = 1:soilvar.nsoi-1
   tk_plus_onehalf(i) = soilvar.tk(i) * soilvar.tk(i+1) * (soilvar.z(i)-soilvar.z(i+1)) / ...
   (soilvar.tk(i)*(soilvar.z_plus_onehalf(i)-soilvar.z(i+1)) + soilvar.tk(i+1)*(soilvar.z(i)-soilvar.z_plus_onehalf(i)));
end

% --- Heat flux at time n (W/m2)

switch solution
   case 'Crank-Nicolson'
   for i = 1:soilvar.nsoi-1
      f(i) = -tk_plus_onehalf(i) * (soilvar.tsoi(i) - soilvar.tsoi(i+1)) / soilvar.dz_plus_onehalf(i);
   end
   f(soilvar.nsoi) = 0;
end

% --- Set up tridiagonal matrix

% Top soil layer with tsurf as boundary condition

i = 1;
m = soilvar.cv(i) * soilvar.dz(i) / dt;
a(i) = 0;
c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
b(i) = m - c(i) + soilvar.tk(i) / (0 - soilvar.z(i));
d(i) = m * soilvar.tsoi(i) + soilvar.tk(i) / (0 - soilvar.z(i)) * tsurf;

switch solution
   case 'Crank-Nicolson'
   a(i) = 0;
   c(i) = -0.5 * tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
   b(i) = m - c(i) + soilvar.tk(i) / (0 - soilvar.z(i));
   d(i) = m * soilvar.tsoi(i) + 0.5 * f(i) + soilvar.tk(i) / (0 - soilvar.z(i)) * tsurf;
end

% Layers 2 to nsoi-1

for i = 2:soilvar.nsoi-1
   m = soilvar.cv(i) * soilvar.dz(i) / dt;
   a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
   c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
   b(i) = m - a(i) - c(i);
   d(i) = m * soilvar.tsoi(i);
end

switch solution
   case 'Crank-Nicolson'
   for i = 2:soilvar.nsoi-1
      m = soilvar.cv(i) * soilvar.dz(i) / dt;
      a(i) = -0.5 * tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
      c(i) = -0.5 * tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
      b(i) = m - a(i) - c(i);
      d(i) = m * soilvar.tsoi(i) + 0.5 * (f(i) - f(i-1));
   end
end

% Bottom soil layer with zero heat flux

i = soilvar.nsoi;
m = soilvar.cv(i) * soilvar.dz(i) / dt;
a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
c(i) = 0;
b(i) = m - a(i);
d(i) = m * soilvar.tsoi(i);

switch solution
   case 'Crank-Nicolson'
   a(i) = -0.5 * tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
   c(i) = 0;
   b(i) = m - a(i);
   d(i) = m * soilvar.tsoi(i) - 0.5 * f(i-1);
end

% --- Solve for soil temperature

[soilvar.tsoi] = tridiagonal_solver (a, b, c, d, soilvar.nsoi);

% --- Derive energy flux into soil (W/m2)

soilvar.gsoi = soilvar.tk(1) * (tsurf - soilvar.tsoi(1)) / (0 - soilvar.z(1));

% --- Phase change for soil layers undergoing freezing of thawing

switch soilvar.method

   case 'apparent-heat-capacity'

   % No explicit phase change energy flux. This is included in the heat capacity.

   soilvar.hfsoi = 0;

   case 'excess-heat'

   % Adjust temperatures for phase change. Freeze or melt ice using energy
   % excess or deficit needed to change temperature to the freezing point.
   % The variable hfsoi is returned as the energy flux from phase change (W/m2).

   [soilvar] = phase_change (physcon, soilvar, dt);

end

% --- Check for energy conservation

% Sum change in energy (W/m2)

edif = 0;
for i = 1:soilvar.nsoi
   edif = edif + soilvar.cv(i) * soilvar.dz(i) * (soilvar.tsoi(i) - tsoi0(i)) / dt;
end

% Error check

err = edif - soilvar.gsoi - soilvar.hfsoi;
if (abs(err) > 1e-03)
   error ('Soil temperature energy conservation error')
end

function [soilvar] = soil_thermal_properties (physcon, soilvar)

% Calculate soil thermal conductivity and heat capacity

% ------------------------------------------------------
% Input
%   physcon.hfus             ! Heat of fusion for water at 0 C (J/kg)
%   physcon.tfrz             ! Freezing point of water (K)
%   physcon.rhowat           ! Density of water (kg/m3)
%   physcon.rhoice           ! Density of ice (kg/m3)
%   soilvar.method           ! Use excess heat or apparent heat capacity for phase change
%   soilvar.nsoi             ! Number of soil layers
%   soilvar.dz               ! Soil layer thickness (m)
%   soilvar.tsoi             ! Soil temperature (K)
%   soilvar.h2osoi_liq       ! Unfrozen water, liquid (kg H2O/m2)
%   soilvar.h2osoi_ice       ! Frozen water, ice (kg H2O/m2)
%
% Input/output
%   soilvar.tk               ! Thermal conducitivty (W/m/K)
%   soilvar.cv               ! Volumetric heat capacity (J/m3/K)
% ------------------------------------------------------

% Temperature range for freezing and thawing (K)

tinc = 0.5;

% Unfrozen and frozen thermal conductivity (W/m/K)

tku = 1.860;
tkf = 2.324;

% Unfrozen and frozen heat capacity (J/m3/K)

cvu = 2.862e06;
cvf = 1.966e06;

for i = 1:soilvar.nsoi

   % --- Volumetric soil water and ice

   watliq = soilvar.h2osoi_liq(i) / (physcon.rhowat * soilvar.dz(i));
   watice = soilvar.h2osoi_ice(i) / (physcon.rhoice * soilvar.dz(i));

   % Heat of fusion (J/m3) - This is equivalent to ql = hfus * (h2osoi_liq + h2osoi_ice) / dz

   ql = physcon.hfus * (physcon.rhowat * watliq + physcon.rhoice * watice);

   % Heat capacity and thermal conductivity

   if (soilvar.tsoi(i) > physcon.tfrz+tinc)
      soilvar.cv(i) = cvu;
      soilvar.tk(i) = tku;
   end

   if (soilvar.tsoi(i) >= physcon.tfrz-tinc & soilvar.tsoi(i) <= physcon.tfrz+tinc)
      switch soilvar.method
         case 'apparent-heat-capacity'
         soilvar.cv(i) = (cvf + cvu) / 2 + ql / (2 * tinc);
         case 'excess heat'
         soilvar.cv(i) = (cvf + cvu) / 2;
      end
      soilvar.tk(i) = tkf + (tku - tkf) * (soilvar.tsoi(i) - physcon.tfrz + tinc) / (2 * tinc);
   end

   if (soilvar.tsoi(i) < physcon.tfrz-tinc)
      soilvar.cv(i) = cvf;
      soilvar.tk(i) = tkf;
   end

end

function [u] = tridiagonal_solver (a, b, c, d, n)

% Solve for U given the set of equations R * U = D, where U is a vector
% of length N, D is a vector of length N, and R is an N x N tridiagonal
% matrix defined by the vectors A, B, C each of length N. A(1) and
% C(N) are undefined and are not referenced.
%
%     |B(1) C(1) ...  ...  ...                     |
%     |A(2) B(2) C(2) ...  ...                     |
% R = |     A(3) B(3) C(3) ...                     |
%     |                    ... A(N-1) B(N-1) C(N-1)|
%     |                    ... ...    A(N)   B(N)  |
%
% The system of equations is written as:
%
%    A_i * U_i-1 + B_i * U_i + C_i * U_i+1 = D_i
%
% for i = 1 to N. The solution is found by rewriting the
% equations so that:
%
%    U_i = F_i - E_i * U_i+1

% --- Forward sweep (1 -> N) to get E and F

e(1) = c(1) / b(1);

for i = 2: 1: n-1
   e(i) = c(i) / (b(i) - a(i) * e(i-1));
end

f(1) = d(1) / b(1);

for i = 2: 1: n
   f(i) = (d(i) - a(i) * f(i-1)) / (b(i) - a(i) * e(i-1));
end

% --- Backward substitution (N -> 1) to solve for U

u(n) = f(n);

for i = n-1: -1: 1
   u(i) = f(i) - e(i) * u(i+1);
end