WSL/SLF GitLab Repository

Commit 70138927 authored by Aschauer's avatar Aschauer
Browse files

pooling rho settling parameters

only overburden of the layer is changing settling resistance
parent fbb81d98
......@@ -52,49 +52,19 @@ def _R_reduction_overburden(
return R_r
@njit
def _R_reduction_density(
rho_layers,
rho_minsettling,
rho_maxsettling,
):
R_r = np.zeros_like(rho_layers)
for i, r in enumerate(rho_layers):
if r <= rho_maxsettling:
R_r[i] = 1
elif r >= rho_minsettling:
R_r[i] = 0
else:
R_r[i] = (r - rho_minsettling)/(rho_maxsettling-rho_minsettling)
return R_r
@njit
def _get_settling_resistance(
rho_layers,
overburden_layers,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
ob_minsettling,
ob_maxsettling,
rho_minsettling,
rho_maxsettling,
):
# calculate weights (scaling factors) of density and overburden.
w_rho = (R_max-R_min) * ratio_settling_influence_rho_vs_ob/(ratio_settling_influence_rho_vs_ob+1)
w_sigma = (R_max-R_min) * 1/(ratio_settling_influence_rho_vs_ob+1)
# calculate reduction of R due to density and overburden.
R_reduction_rho = _R_reduction_density(rho_layers, rho_minsettling, rho_maxsettling)
R_reduction_sigma = _R_reduction_overburden(overburden_layers, ob_minsettling, ob_maxsettling)
R = R_max - w_rho*R_reduction_rho - w_sigma*R_reduction_sigma
R_reduction_sigma = _R_reduction_overburden(
overburden_layers, ob_minsettling, ob_maxsettling)
R = R_max - (R_max-R_min)*R_reduction_sigma
return R
......@@ -109,7 +79,7 @@ def _adjust_rho_max_wet_snowpack(
swe_layers > 0,
rho_max_wet - (rho_max_wet-rho_max_layers)*np.exp(-wetting_speed),
rho_max_layers
)
)
@njit
......@@ -146,20 +116,21 @@ def _adjust_rho_max_based_on_overburden(
:class:`numpy.ndarray`
Adapted rho_max in the layers based on overburden
"""
rho_max_current = np.zeros_like(overburden_layers)
for i, o in enumerate(overburden_layers):
if o >= max_overburden:
rho_max_current[i] = rho_max_wet
else:
rho_max_current[i] = ((rho_max_wet-rho_max_dry)/max_overburden)*o + rho_max_dry
updated_rho_max_layers = np.where((rho_max_current>rho_max_layers),
rho_max_current,
rho_max_current[i] = ((rho_max_wet-rho_max_dry) /
max_overburden)*o + rho_max_dry
updated_rho_max_layers = np.where((rho_max_current > rho_max_layers),
rho_max_current,
rho_max_layers)
updated_rho_max_layers = np.where((swe_layers>0), updated_rho_max_layers, 0.)
updated_rho_max_layers = np.where((swe_layers > 0), updated_rho_max_layers, 0.)
return updated_rho_max_layers
......@@ -185,11 +156,11 @@ def _remove_swe_from_top(
swe_removed = 0
l = len(swe_layers)-1
# melting from top means going backward in layers.
while swe_removed > delta_swe: # both are (or will be) negative
while swe_removed > delta_swe: # both are (or will be) negative
swe_removed = swe_removed - swe_layers[l]
swe_layers[l] = 0
l = l-1
# minimal floating point errors can cause the while loop to
# minimal floating point errors can cause the while loop to
# run away, in that case we force it to stop at bottom of the
# snowpack.
if l == -1:
......@@ -212,13 +183,10 @@ def timestep_forward(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
wetting_speed,
wetting_speed,
):
"""Process a timestep forward.
......@@ -287,27 +255,24 @@ def timestep_forward(
rho_max_dry,
rho_max_wet,
max_overburden
)
)
# calculate settling resistance
R = _get_settling_resistance(
rho_layers,
overburden_layers,
R_max=R_max,
R_min=R_min,
ratio_settling_influence_rho_vs_ob=ratio_settling_influence_rho_vs_ob,
ob_minsettling=overburden_minsettling,
ob_maxsettling=overburden_maxsettling,
rho_minsettling=rho_minsettling,
rho_maxsettling=rho_maxsettling,
)
)
# update rho, i.e. calculate settling:
rho_layers = np.where(
(swe_layers>0),
(swe_layers > 0),
rho_max_layers - (rho_max_layers-rho_layers) * np.exp(-1/R),
0.)
if delta_swe > 0:
# rho_new should not get modified in the first timestep.
rho_layers[-1] = rho_new
......@@ -355,9 +320,6 @@ def swe2hs_snowpack_evolution(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
......@@ -365,13 +327,12 @@ def swe2hs_snowpack_evolution(
):
"""
Snowpack evolution within the swe2hs model.
Meant to be called on single hydrological years or chunks of nonzeros
"""
assert R_max >= R_min
assert overburden_minsettling < overburden_maxsettling
assert rho_maxsettling < rho_minsettling
# allocate output arrays:
hs_out = np.zeros(len(swe_input))
......@@ -380,15 +341,15 @@ def swe2hs_snowpack_evolution(
rho_max_layers_out = np.zeros((len(swe_input), len(swe_input)))
# allocate layer containers
swe_layers = np.zeros(0) # tracking of swe
rho_layers = np.zeros(0) # tracking of density
rho_max_layers = np.zeros(0) #tracking of maximum density
swe_layers = np.zeros(0) # tracking of swe
rho_layers = np.zeros(0) # tracking of density
rho_max_layers = np.zeros(0) # tracking of maximum density
# iterate through input array.
for i, swe in enumerate(swe_input):
# get delta swe:
if i==0:
if i == 0:
delta_swe = swe
else:
delta_swe = swe - swe_input[i-1]
......@@ -409,9 +370,6 @@ def swe2hs_snowpack_evolution(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
......@@ -421,9 +379,9 @@ def swe2hs_snowpack_evolution(
hs_layers = _calculate_hs_layers(swe_layers, rho_layers)
hs_out[i] = np.sum(hs_layers)
hs_layers_out[:i+1, i] = hs_layers
rho_max_layers_out[:i+1, i] = rho_max_layers
rho_layers_out[:i+1, i] = rho_layers
......
......@@ -69,9 +69,6 @@ def _swe2hs_on_nonzero_chunks(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
......@@ -90,9 +87,6 @@ def _swe2hs_on_nonzero_chunks(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
......@@ -114,9 +108,6 @@ def swe2hs_1d(
rho_max_wet=RHO_MAX_WET,
R_max=R_MAX,
R_min=R_MIN,
ratio_settling_influence_rho_vs_ob=RATIO_SETTLING_INFLUENCE_RHO_VS_OB,
rho_minsettling=RHO_MINSETTLING,
rho_maxsettling=RHO_MAXSETTLING,
overburden_minsettling=OVERBURDEN_MINSETTLING,
overburden_maxsettling=OVERBURDEN_MAXSETTLING,
max_overburden=MAX_OVERBURDEN,
......@@ -292,9 +283,6 @@ def swe2hs_1d(
rho_max_wet,
R_max,
R_min,
ratio_settling_influence_rho_vs_ob,
rho_minsettling,
rho_maxsettling,
overburden_minsettling,
overburden_maxsettling,
max_overburden,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment