-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathInitialSurfaceMelt_2ColumnFirn.m
138 lines (121 loc) · 6.15 KB
/
InitialSurfaceMelt_2ColumnFirn.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
function [Temp_Profile, rho, k, cp, T_sfcOLD, Water_Content, Sfrac, MeltHours, Lens_Count, S, dHdtREC, dHdtTOT, ksi, TotalMelt, MeltTime, Remaining_pore_space, MeltIn]=InitialSurfaceMelt_2ColumnFirn(t, Temp_Profile, tsm, rho, T_sfcOLD, Water_Content, Sfrac, q, MeltHours, Lens_Count, S, M, dHdtREC, dHdtTOT, ksi, TotalMelt, MeltTime, Remaining_pore_space, Loop, Lens_Level, I, MeltBegins, T_air);
%This function calculates the transfer of heat through the ice and
%determines if surface melting is taking place. If this does occur the
%meltwater is allowed to percolate down through the snow, refreezing once
%it reaches a layer that is below freezing temperature. Ice lens formation
%is monitored.
m=0; % For use in the pde solver. m=0 because the problem is a slab.
Lfus=334000;% Latent heat of fusion
rho_ice=917; %Density of ice
rho_liq=1000; %Density of water
rho_imp=830; %Pore closure depth: firn with the density or above is impermeable to water
k_air=0.022;
cp_air=1004;
k_water=0.5818;
cp_water=4217;
Exposed_Water_Time=0;
MeltIn=0;
for i=1:M %Calculate k and cp for the solid part of each layer
if Temp_Profile(i)>273.15
cp_ice(i)=4186.8;
k_ice(i)=1000*(1.017*10^(-4)+1.695*10^(-6)*Temp_Profile(i));
else
k_ice(i)=1000*(2.24*10^(-3)+5.975*10^(-6)*(273.15-Temp_Profile(i))^1.156); %
cp_ice(i)=1000*(7.16*10^(-3)*Temp_Profile(i)+0.138);%Alexiades & Solomon pg 8
end
end
for i=1:M %Calculate the solid fraction of each layer (from 0 being no ice to 1 being completely ice)
Sfrac(i)=(rho(i)-(Water_Content(i)*rho_ice))/rho_ice; %Needs to be rho_ice for the water part here as it's taken as a fraction of solid ice. Assumption that the densities are close or things get too complicated in the refreezing. Area for possible future development!
if Sfrac(i)<0
Sfrac(i)=rho(i)/rho_ice;
end
end
%Calculate the overall k and cp for each layer
for i=1:M
Air(i)=(1-(Sfrac(i)+Water_Content(i)));
end
k=Sfrac.*k_ice+Air.*k_air+Water_Content'.*k_water;
cp=Sfrac.*cp_ice+Air.*cp_air+Water_Content'.*cp_water;
%Solve the surface energy balance equation for this timestep
Temp_Profile_fit=griddedInterpolant(ksi',Temp_Profile,'pchip');
tbcfun = @(x,q, Temp_Profile_fit, S) -0.98*5.670373*(10^-8)*x^4+q-k(end,end)*(-Temp_Profile_fit(0.999)+x)/(S/1000); % The parameterized function.
T_sfcOLD = fzero(@(x) tbcfun(x,q, Temp_Profile_fit, S), T_air);
%Is there melt?
if T_sfcOLD>=273.15 && q>0
kdTdz=(273.15-Temp_Profile_fit(0.999))*abs(k(end))/(S/1000);
dHdt=(q-0.98*5.670373*(10^-8)*273.15^4-kdTdz)/(Sfrac(end)*Lfus*rho_ice)*tsm;
if dHdt<0 %Error checking for instability
disp('Error- melting at top calculation')
return
end
MeltHours=MeltHours+1;
T_sfcOLD=273.15;
else
dHdt=0;
end
dHdtREC(Loop)=dHdt;%Store for output later
dHdtTOT(Loop+1)=dHdtTOT(Loop)+dHdt;%Records total height lost from surface due to melt
%Recalculate temp and density profiles taking off top melt.
if dHdt>0
[rho, Water_Content]=Profile_After_Melt(rho, dHdt, M, ksi, S, Water_Content);
S=S-dHdt;
%Calculate what happens to the meltwater produced
MeltWaterVol=dHdt*Sfrac(end)*rho_ice/rho_liq;
TotalMelt=TotalMelt+MeltWaterVol;
MeltTime(Loop)=MeltWaterVol;
if Lens_Count>0
MeltIn=MeltWaterVol;
MeltWaterVol=0;%All meltwater is transported into the lake if an ice lens is present
end
if Lens_Count==0 %If no ice lens is present water is allowed to percolate through the firn until it refreezes
[Temp_Profile, rho, I, Water_Content,Remaining_pore_space, Lens_Count, MeltWaterVol, Lens_Level]=Percolation(Temp_Profile, rho, ksi, dHdt, M, S, Sfrac, cp, Water_Content, Lens_Count, Lens_Level, Remaining_pore_space, MeltWaterVol);
end
%Recalculate k and cp for these new profiles
for i=1:M
if Temp_Profile(i)>273.15
cp_ice(i)=4186.8;
k_ice(i)=1000*(1.017*10^(-4)+1.695*10^(-6)*Temp_Profile(i));
else
k_ice(i)=1000*(2.24*10^(-3)+5.975*10^(-6)*(273.15-Temp_Profile(i))^1.156); %
cp_ice(i)=1000*(7.16*10^(-3)*Temp_Profile(i)+0.138);%Alexiades & Solomon pg 8
end
Sfrac(i)=(rho(i)-(Water_Content(i)*rho_ice))/rho_ice; %Needs to be rho_ice for the water part here as it's taken as a fraction of solid ice. Assumption that the densities are close or things get too complicated in the refreezing. Area for possible future development!
if Sfrac(i)<0
Sfrac(i)=rho(i)/rho_ice;
end
end
k=Sfrac.*k_ice+(1-Sfrac).*k_air;%Calculate overall k and cp for each layer
cp=Sfrac.*cp_ice+(1-Sfrac).*cp_air;
end
%Fit gridded interpolants to profiles to allow them to be put into the pde solver
Temp_Profile_fit=griddedInterpolant(ksi',Temp_Profile,'pchip');
rho_fit=griddedInterpolant(ksi',rho,'pchip');
k_fit=griddedInterpolant(ksi',k,'pchip');
cp_fit=griddedInterpolant(ksi',cp,'pchip');
%Use the pde solver pdepe to calculate heat transfer
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,ksi,t,[],rho_fit, cp_fit, k_fit, S, Temp_Profile_fit, q, T_sfcOLD, dHdt);
Temperature = sol(:,:,1);%Extract the first solution component;
Temp_Profile=Temperature(end,:).';
for i=1:M %Error checking to check for stability of pde solver
if Temp_Profile(i)~=real( Temp_Profile(i)) || Temp_Profile(i)>300 || Temp_Profile(i)<200
disp ('Error pde solution unstable')
break
end
end
%%Pde heat transfer calculations pre-exposed water
function [c,f,s] = pdex1pde(ksi,t,u,DuDx,rho_fit, cp_fit, k_fit, S, Temp_Profile_fit, Q, T_sfcOLD, dHdt)
%global rho_fit cp_fit k_fit S
c = cp_fit(ksi)*rho_fit(ksi)*(S^2);
f = k_fit(ksi)*DuDx;
s = S*ksi*DuDx*(dHdt/3600)*cp_fit(ksi)*rho_fit(ksi);
% --------------------------------------------------------------
function u = pdex1ic(ksi,rho_fit, cp_fit, k_fit, S, Temp_Profile_fit, Q, T_sfcOLD, dHdt)
%global Temp_Profile_fit
u = Temp_Profile_fit(ksi);
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t,rho_fit, cp_fit, k_fit, S, Temp_Profile_fit, Q, T_sfcOLD, dHdt)
%global Q T_sfcOLD
pr = ur-T_sfcOLD;%Sets surface temperature to that calculated by solving the surface energy balance equation
qr = 0;
pl = 0;
ql = 1; %sets right side temperature grad to 0