-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathS3M_2D_physics.py
229 lines (181 loc) · 10.3 KB
/
S3M_2D_physics.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
# Re-writing of a punctual version of S3M with python language. Project related to the doctoral thesis:Data assimilation
# and deep learning·I will try to re-write S3M in python and develop an EnKF or PF data assimilation procedure to
# replace the current nudging procedure. Then I will use data assimilation data to train a neural network (CNN/LSTM/RNN)
# to exploit the computational capacity of DL to make this procedure comparable with operative times.
# In this version I will omit the glacial components.
# #----------------------------------------------------------------------------------------------------------------------
# ## Input data of S3M are:
# 1) Land data, (mandatory but not for punctual model)
# 2) Meteorological observations, (mandatory)
# 3) SCA satellite images, (optional)
# 4) SWE independent estimates for assimilation, (optional)
# 5) Ice thickness and a number of other ancillary glacier data, (optional).[ This will not be needed]
# ## Output results are (among others):
# 1) Snowpack runoff,
# 2) SWE (dry and wet), snow density (dry and wet), snow depth, bulk liquid water content
# 3) Snowfall, Rainfall, and Precipitation rates, as well as fresh-snow density
# 4) Snow age, albedo, melt, and snowpack runoff
# 5) Ice thickness. [ This will not be needed]
# -------------------------------------------------------------------------------------
# -------------------------------------------------------------------------------------
# Library
from lib_utilis_flux import PhasePart, density, melting, refreezing, Hydraulics
from solar_radiation import solar_radiation, solarhours
import numpy as np
def S3M_2D_physics(log_stream, meteo, dynamic_inputs, parameters, state_vector, output_vector, Time, change_part):
# now every variable has to be a vector 40x1 (40 measurements points)
mass_balance = 0
Outflow_ExcessRain = np.zeros(40) # Initialize as a NumPy array
Outflow_ExcessMelt = np.zeros(40) # Initialize as a NumPy array
# Meteorological input upload
T_air = meteo[:,dynamic_inputs["T_air"]]
P = meteo[:, dynamic_inputs["precipitation"]]
RH = meteo[:,dynamic_inputs["Relative_Humidity"]]
lat = parameters["latitude"] * np.ones(40)
lon = parameters["longitude"] * np.ones(40)
# Day of the year and hour of the day from the timestamp,but Time is a vector of 40 elements, so also h and doy
h = Time.hour
doy = Time.timetuple().tm_yday
# apply the solar radiation function to the matrix of latitudes and longitudes
Rtoa = solar_radiation(lat, lon, doy, h) # Rtoa is a vector of 40 elements
hrise, hset = solarhours(lat, lon, doy)
# Calculate the solar radiation for the day , if the hour is between sunrise and sunset the min between the actual
# solar radiation and the maximum solar radiation is taken, otherwise 0. consider thar Rtoa is a vector of 40 elements
#and so is the output vector and hrise and hset
Radiation = np.minimum(Rtoa, meteo[:, dynamic_inputs["Radiation"]])
mask = (h < hrise) | (h > hset)
Radiation[mask] = 0
# ----------------------------------------------------------------------------------------
# first sanity check
# Vectorized operation to set output_vector[10] to 0 if it is less than 0.01
mask = output_vector[:, 10] < 0.01
output_vector[mask, 10] = 0
# Vectorized operation to set state_vector elements to 0
state_vector[mask, :-1] = 0
# Vectorized operation to set output_vector elements to 0 from index 11 onwards
output_vector[mask, 11:] = 0
SWE_W, SWE_D = state_vector[:, 0], state_vector[:,1]
Sf_daily_cum = output_vector[:, 5]
# ------------------------------------------------------------------------------------
# Precipitation phase partitioning
alpha, beta, gamma = parameters["alpha"], parameters["beta"], parameters["gamma"]
Snowfall, Rainfall = PhasePart(P, alpha, beta, gamma, T_air, RH, change_part)
Sf_daily_cum += Snowfall
# ---------------------------------------------------------------------------------------------
# Vectorized operation to handle Rainfall and SWE_D
mask_rainfall = Rainfall > 0
mask_swe_d = SWE_D >= 10.0
# Update SWE_W where Rainfall > 0 and SWE_D >= 10.0
SWE_W[mask_rainfall & mask_swe_d] += Rainfall[mask_rainfall & mask_swe_d]
# Update Outflow_ExcessRain and SWE_W where Rainfall > 0 and SWE_D < 10.0
mask_swe_d_less = ~mask_swe_d
Outflow_ExcessRain[mask_rainfall & mask_swe_d_less] = Rainfall[mask_rainfall & mask_swe_d_less] + SWE_W[
mask_rainfall & mask_swe_d_less]
SWE_W[mask_rainfall & mask_swe_d_less] = 0.0
# Ensure SWE_D and SWE_W are non-negative
SWE_D = np.maximum(SWE_D, 0)
SWE_W = np.maximum(SWE_W, 0)
# Calculate SWE
SWE = SWE_D + SWE_W
# Compute snow density
Rho_D_min, Rho_D_max, Rho_S_max, RhoW, dt = parameters["RhoSnowMin"], parameters["RhoSnowMax"], parameters[
"RhoFreshSnowMax"], parameters["RhoW"], parameters["dt"]
Rho_D, RhoS0, SnowTemp, H_D = density(Rho_D_min, Rho_D_max, Rho_S_max, RhoW, dt, state_vector, output_vector, SWE_D,
Snowfall, T_air)
# Compute melting and refreezing
cm = dt / 86400
T_10D, T_1D, Ttau, mrad0, mr0 = meteo[:,4], meteo[:,5], parameters["Ttau"], parameters["mrad0"], parameters["mr0"]
As, albedo, multiplicative_term = output_vector[:,11], state_vector[:,3], parameters["multiplicative_albedo"]
Melting, albedo, As, Sf_daily_cum, mrad, mr = melting(Time, mrad0, mr0, T_air, T_10D, T_1D, Ttau, Radiation, RhoW,
dt, cm, SWE_D, albedo, As, SWE, Sf_daily_cum,
multiplicative_term)
Refreezing = refreezing(T_air, T_10D, SWE_W, mr0, cm, Ttau)
# -------------------------------------------------------------------------------------
# Ensure Melting and Refreezing are non-negative
Melting = np.maximum(Melting, 0)
Refreezing = np.maximum(Refreezing, 0)
# Vectorized operation for Melting
mask_melting = (0 < SWE_D) & (SWE_D <= Melting)
Melting[mask_melting] = SWE_D[mask_melting]
Outflow_ExcessMelt[mask_melting] = SWE_D[mask_melting] + SWE_W[mask_melting]
SWE_D[mask_melting] = 0.0
SWE_W[mask_melting] = 0.0
SWE[mask_melting] = 0.0
mask_no_melting = ~mask_melting
SWE_D[mask_no_melting] -= Melting[mask_no_melting]
SWE_W[mask_no_melting] += Melting[mask_no_melting]
SWE[mask_no_melting] = SWE_D[mask_no_melting] + SWE_W[mask_no_melting]
# Vectorized operation for Refreezing
mask_refreezing = (0 < SWE_W) & (SWE_W <= Refreezing)
Rho_D[mask_refreezing] = (SWE_D[mask_refreezing] + Refreezing[mask_refreezing]) / \
((Refreezing[mask_refreezing] / 917) + (SWE_D[mask_refreezing] / Rho_D[mask_refreezing]))
SWE_D[mask_refreezing] += Refreezing[mask_refreezing]
SWE_W[mask_refreezing] = 0.0
SWE[mask_refreezing] = SWE_D[mask_refreezing]
mask_partial_refreezing = (0 < Refreezing) & (Refreezing < SWE_W)
Rho_D[mask_partial_refreezing] = (SWE_D[mask_partial_refreezing] + Refreezing[mask_partial_refreezing]) / \
((Refreezing[mask_partial_refreezing] / 917) + (
SWE_D[mask_partial_refreezing] / Rho_D[mask_partial_refreezing]))
SWE_D[mask_partial_refreezing] += Refreezing[mask_partial_refreezing]
SWE_W[mask_partial_refreezing] -= Refreezing[mask_partial_refreezing]
SWE[mask_partial_refreezing] = SWE_D[mask_partial_refreezing] + SWE_W[mask_partial_refreezing]
# Ensure SWE_D and SWE_W are non-negative
SWE_D = np.maximum(SWE_D, 0)
SWE_W = np.maximum(SWE_W, 0)
SWE = SWE_D + SWE_W
# Update height
H_D = ((SWE_D / 1000) * RhoW) / Rho_D
# Outflow
outflow, H_S = Hydraulics(Rho_D, RhoW, SWE_D, SWE_W, H_D, dt)
# Ensure SWE_W is non-negative and update SWE
mask_swe_w_positive = SWE_W > 0
SWE_W[mask_swe_w_positive] -= outflow[mask_swe_w_positive]
SWE[mask_swe_w_positive] = SWE_D[mask_swe_w_positive] + SWE_W[mask_swe_w_positive]
# Update outflow
outflow += Outflow_ExcessMelt + Outflow_ExcessRain
# Handle cases where SWE_W is less than or equal to 0
mask_swe_w_non_positive = SWE_W <= 0
outflow[mask_swe_w_non_positive] += SWE_W[mask_swe_w_non_positive]
SWE_W[mask_swe_w_non_positive] = 0.0
SWE[mask_swe_w_non_positive] = SWE_D[mask_swe_w_non_positive]
# Log error if SWE_W is negative
mask_swe_w_negative = SWE_W < 0
if np.any(mask_swe_w_negative):
log_stream.error('negative swe_w ' + str(SWE_W[mask_swe_w_negative]) + ' at time ' + str(Time))
# Handle cases where SWE is less than 0
mask_swe_negative = SWE < 0
if np.any(mask_swe_negative):
log_stream.error('negative swe ' + str(SWE[mask_swe_negative]) + ' at time ' + str(Time))
SWE[mask_swe_negative] = 0
# Update height
H_D = ((SWE_D / 1000) * RhoW) / Rho_D
H_S = np.where((SWE_W / 1000) >= ((1 - Rho_D / 917) * H_D),
H_D + (SWE_W / 1000) - (1 - Rho_D / 917) * H_D,
H_D)
# Calculate theta_w and Rho_s
mask_h_s_positive = H_S > 0
theta_w = np.where(mask_h_s_positive, SWE_W / 1000 / H_S, 0)
Rho_s = np.where(mask_h_s_positive, (Rho_D * H_D + RhoW * (SWE_W / 1000)) / H_S, 0)
# Mass Balance
mass_balance = np.where(
(np.round((SWE - output_vector[:, 10]), 2) != np.round((Snowfall + Rainfall - outflow), 2)) &
(np.round((SWE - output_vector[:, 10]), 1) != np.round((Snowfall + Rainfall - outflow), 1)),
1, 0
)
# Update Time
Time = Time.timestamp()
# Create new vectors
input_vector_new = meteo
variables = [SWE_W, SWE_D, Rho_D, albedo]
state_vector_new = np.zeros((40, 4))
# Loop through the variables and assign them to state_vector
for idx, var in enumerate(variables):
state_vector_new [:, idx] = var
# Create new output vector
output_variables=[Rainfall, Snowfall, Melting, Refreezing, outflow, Sf_daily_cum, Time, mass_balance,
mrad, mr, SWE, As, H_D, theta_w, H_S, Rho_s]
output_vector_new = np.zeros((40, 16))
# Loop through the output variables and assign them to output_vector
for idx, var in enumerate(output_variables):
output_vector_new[:, idx] = var
return meteo, state_vector_new, output_vector_new, mass_balance