-
Notifications
You must be signed in to change notification settings - Fork 0
/
msm_old.py
119 lines (97 loc) · 4.52 KB
/
msm_old.py
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
"""
Multilayer Subsurface Model (of a glacier).
Simulates surface temperature changes according to supra-glacial heat balance.
Yields in-glacier heat flux and a flux, needed to warp the surface up to zero degrees Celsius.
Based on equations from:
(Lisette) Klok, E., & Oerlemans, J. (2002). Model study of the spatial distribution of the energy
and mass balance of Morteratschgletscher, Switzerland. Journal of Glaciology, 48(163), 505-518.
doi:10.3189/172756502781831133
"""
from dataclasses import dataclass, field
from typing import List
from parameter_classes import CONST
import csv
import matplotlib.pyplot as plt
@dataclass
class Layer:
depth: float # meters
temp: float # degree Celsius
@dataclass
class Msm:
thickness_list: List[float]
temp_list: List[float]
layers: List[Layer] = field(init=False)
c: float = field(init=False) # specific heat capacity of ice
k: float = field(init=False) # thermal diffusivity
density: float = field(init=False)
def __post_init__(self):
# initializing constants:
self.c = CONST["specific_heat_capacity_ice"]
self.k = CONST["thermal_diffusivity_ice"]
self.density = CONST["ice_density"]
# initializing layers:
self.layers = []
if len(self.thickness_list) and len(self.temp_list):
for i in range(0, len(self.thickness_list)): # TODO: use zip() here
self.layers.append(Layer(self.thickness_list[i], self.temp_list[i]))
def tick(self, timestep, atmo_heat_flux, verbose=False):
grad_t = (self.layers[1].temp - self.layers[0].temp) / (self.layers[1].depth - self.layers[0].depth)
print(grad_t)
glacier_heat_flux = self.k * grad_t * self.c * self.density
delta_temp_surf = self.k * grad_t / self.layers[1].depth + atmo_heat_flux / (self.c * self.density * self.layers[1].depth)
grad_t_subs = (self.layers[2].temp - self.layers[1].temp) / (self.layers[2].depth - self.layers[1].depth)
print(grad_t_subs)
delta_temp_subs = self.k * (grad_t - grad_t_subs) / self.layers[2].depth
# print(delta_temp_subs)
full_balance = atmo_heat_flux + glacier_heat_flux
# flux, needed to warm the surface up to 0 degrees:
zero_flux = -self.layers[0].temp * self.c * self.density * self.layers[1].depth / timestep
# melt_flux = 0 if full_balance <= zero_flux else full_balance - zero_flux
# print("Melt flux is %.1f W/m^2" % melt_flux)
# return melt_flux
self.layers[0].temp += delta_temp_surf * timestep
self.layers[0].temp = 0 if self.layers[0].temp > 0 else self.layers[0].temp # temperature can't exceed melting temp
self.layers[1].temp += delta_temp_subs * timestep
self.layers[1].temp = 0 if self.layers[1].temp > 0 else self.layers[1].temp
# the last, most deep layer, has an unlimited heat capacity, so its temperature cannot be changed
if verbose:
print("Msm timestep is %s seconds" % timestep)
print("Glacier heat flux is %.1f W/m^2" % glacier_heat_flux)
print("Temperature change rate is %.6f/s or %.1f/step" % (delta_temp_surf, delta_temp_surf * timestep))
print("Full heat balance is %.1f W/m^2" % full_balance)
print("Zero-flux is %.1f W/m^2" % zero_flux)
return glacier_heat_flux, zero_flux
if __name__ == "__main__":
m = Msm([0, 0.22, 2.78], [0, -0.15, -2.5])
print(m)
t_surf = []
t_2 = []
t_3 = []
# balances = []
ground_heat_fluxes = []
# with open("/home/tepex/PycharmProjects/energy/aws/out_20190806_full.csv") as csvfile:
with open("/home/tepex/PycharmProjects/energy/aws/out_Aldegonda_hourly_data_GMT_total.csv") as csvfile:
reader = csv.DictReader(csvfile)
for row in reader:
# print(row["TOTAL"])
# g, zero_f = m.tick(3600, float(row["TOTAL_BALANCE"]))
g, zero_f = m.tick(3600, 0)
# print(m)
# balances.append(float(row["TOTAL_BALANCE"]))
ground_heat_fluxes.append(g)
t_surf.append(m.layers[0].temp)
t_2.append(m.layers[1].temp)
t_3.append(m.layers[2].temp)
plt.style.use("seaborn")
plt.figure(figsize=(5, 3))
plt.plot(t_surf, label="glacier surface")
plt.plot(t_2, label="subsurface layer")
plt.plot(t_3, label="glacier body")
plt.legend()
plt.tight_layout()
plt.show()
plt.figure(figsize=(5, 3))
plt.plot(ground_heat_fluxes, label="ground heat flux")
plt.legend()
plt.tight_layout()
plt.show()