-
Notifications
You must be signed in to change notification settings - Fork 0
/
model.py
152 lines (119 loc) · 5.12 KB
/
model.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
import numpy as np
import matplotlib.pylab as plt
import math
from utils import read_data
class VariablePrep(object):
def __init__(self, filename="LATTERBACHmeteo.txt", cal_factor=0.8, cal_data=True, val_data=False):
self.area = 562 # in km2
self.elevation_basin = 1590 # mean elevation of the catchment
self.elevation_station = 1320 # elevation of the measuring station
self.lapse_rate = 0.5/100 # lapse rate for temperature in Celsius/metre
assert val_data != cal_data, "Select one of cal_data or val_data"
j_day, precip, t_max, t_min, q_obs = read_data(filename)
div_index = int(np.floor(cal_factor*j_day.shape[0]))
if cal_data:
self.j_day = j_day[:div_index]
self.precip = precip[:div_index]
self.t_max = t_max[:div_index]
self.t_min = t_min[:div_index]
self.q_obs = q_obs[:div_index]
elif val_data:
self.j_day = j_day[div_index:]
self.precip = precip[div_index:]
self.t_max = t_max[div_index:]
self.t_min = t_min[div_index:]
self.q_obs = q_obs[div_index:]
self.q = self.q_obs*60*60*24*1000/(self.area*1e6)
self.t_min_av = (self.elevation_station -
self.elevation_basin)*self.lapse_rate + self.t_min
self.t_max_av = (self.elevation_station -
self.elevation_basin)*self.lapse_rate + self.t_max
self.t_av = (self.t_min_av+self.t_max_av)/2
class ConceptualWatbalModel(object):
def __init__(self, parameters):
self.snow_cover_prev = 0
self.ss_prev = 10
self.sg_prev = 100
self.k = parameters["k"]
self.s_max = parameters["s_max"]
self.fr = parameters["fr"]
self.rg = parameters["rg"]
def divide_snow_rain(self):
if self.temp_max <= 0:
self.snow = self.precip
self.rainfall = 0
if self.temp_min > 0:
self.snow = 0
self.rainfall = self.precip
if (self.temp_max > 0 and self.temp_min <= 0):
self.rainfall = self.temp_max / \
(self.temp_max-self.temp_min)*self.precip
self.snow = self.precip - self.rainfall
def compute_snow_melt(self, temp_snow_melt=0):
self.temp_snow_melt = temp_snow_melt
if self.temp <= self.temp_snow_melt:
self.melt = 0
elif self.temp > self.temp_snow_melt:
self.melt = min(
self.k*( self.temp-self.temp_snow_melt), self.snow_cover_prev)
def compute_snow_cover(self):
self.snow_cover_current = self.snow_cover_prev - self.melt + self.snow
def compute_pot_et(self, lat=46.5):
phi = lat
delta = 0.4093*np.sin((2*math.pi/365)*self.j_day-1.405)
omega_s = np.arccos(-np.tan(2*np.pi*phi/360)*np.tan(delta))
Nt = 24*omega_s/np.pi
a, b, c = 0.6108, 17.27, 237.3
es = a*np.exp(b*self.temp/(self.temp+c))
E = (2.1*(Nt**2)*es/(self.temp+273.3))
if self.temp <= 0:
E = 0
self.pot_et = E
def compute_et(self):
self.et = self.ss_prev/self.s_max*self.pot_et
def compute_soil_zone_water_content(self):
self.ss_current = self.ss_prev + self.melt+self.rainfall-self.et
if self.ss_current < 0:
self.ss_current = 0
def compute_surface_runoff(self):
if self.ss_current > self.s_max:
self.q_surf = (self.rainfall+self.melt)*self.fr
self.ss_current = self.ss_current - self.q_surf
else:
self.q_surf = 0
def compute_percolation_groundwater(self):
if self.ss_current > self.s_max:
self.perc_gw = self.ss_current - self.s_max
self.ss_current = self.ss_current-self.perc_gw
else:
self.perc_gw = 0
def compute_ground_water_reservoir(self):
self.q_gw = (1/self.rg)*self.sg_prev
self.sg_current = self.sg_prev + self.perc_gw - self.q_gw
if self.sg_current < 0:
self.q_gw = 0
self.sg_current = 0
def compute_stream_flow(self):
self.q_out = self.q_surf + self.q_gw
def update(self):
self.snow_cover_prev = self.snow_cover_current
self.ss_prev = self.ss_current
self.sg_prev = self.sg_current
def simulate(self, variables):
self.j_day = variables["j_day"]
self.temp_max = variables["temp_max"]
self.temp_min = variables["temp_min"]
self.temp = variables["temp"]
self.precip = variables["precip"]
self.divide_snow_rain()
self.compute_snow_melt()
self.compute_snow_cover()
self.compute_pot_et()
self.compute_et()
self.compute_soil_zone_water_content()
self.compute_surface_runoff()
self.compute_percolation_groundwater()
self.compute_ground_water_reservoir()
self.compute_stream_flow()
self.update()
return self.q_out, self.et, self.pot_et, self.precip, self.snow, self.rainfall, self.melt, self.snow_cover_current, self.ss_current, self.q_surf, self.sg_current, self.q_gw, self.temp_min, self.temp, self.temp_max