GR4J Model¶
GR4J (Genie Rural a 4 parametres Journalier) is a daily lumped rainfall-runoff model with four free parameters. It combines a production store, a routing store, and two unit hydrographs (UH1/UH2) for routing.
- Registered name:
gr4j(seeMODEL_DICT) - Parameters: 4
- Routing: unit hydrographs + nonlinear routing store
API Reference¶
calculate_evap_store(s, evap_net, x1)
¶
Determines the evaporation loss from the production store
Source code in hydromodel/models/gr4j.py
@jit(nopython=True)
def calculate_evap_store(s, evap_net, x1):
"""Determines the evaporation loss from the production store"""
n = s * (2.0 - s / x1) * np.tanh(evap_net / x1)
d = 1.0 + (1.0 - s / x1) * np.tanh(evap_net / x1)
return n / d
calculate_perc(current_store, x1)
¶
Determines how much water percolates out of the production store to streamflow
Source code in hydromodel/models/gr4j.py
@jit(nopython=True)
def calculate_perc(current_store, x1):
"""Determines how much water percolates out of the production store to streamflow"""
return current_store * (
1.0 - (1.0 + (4.0 / 9.0 * current_store / x1) ** 4) ** -0.25
)
calculate_precip_store(s, precip_net, x1)
¶
Calculates the amount of rainfall which enters the storage reservoir.
Source code in hydromodel/models/gr4j.py
@jit(nopython=True)
def calculate_precip_store(s, precip_net, x1):
"""Calculates the amount of rainfall which enters the storage reservoir."""
n = x1 * (1.0 - (s / x1) ** 2) * np.tanh(precip_net / x1)
d = 1.0 + (s / x1) * np.tanh(precip_net / x1)
return n / d
gr4j(p_and_e, parameters, warmup_length, return_state=False, normalized_params='auto', **kwargs)
¶
run GR4J model
Parameters¶
ndarray
3-dim input -- [time, basin, variable]: precipitation and potential evaporation
parameters 2-dim variable -- [basin, parameter]: the parameters are x1, x2, x3 and x4 warmup_length length of warmup period return_state if True, return state values, mainly for warmup periods normalized_params parameter format specification: - "auto": automatically detect if parameters are normalized (0-1) or original scale (default) - True: parameters are normalized (0-1 range), will be converted to original scale - False: parameters are already in original scale, use as-is
Returns¶
Union[np.array, tuple] streamflow or (streamflow, states)
Source code in hydromodel/models/gr4j.py
def gr4j(
p_and_e,
parameters,
warmup_length: int,
return_state=False,
normalized_params="auto",
**kwargs,
):
"""
run GR4J model
Parameters
----------
p_and_e: ndarray
3-dim input -- [time, basin, variable]: precipitation and potential evaporation
parameters
2-dim variable -- [basin, parameter]:
the parameters are x1, x2, x3 and x4
warmup_length
length of warmup period
return_state
if True, return state values, mainly for warmup periods
normalized_params
parameter format specification:
- "auto": automatically detect if parameters are normalized (0-1) or original scale (default)
- True: parameters are normalized (0-1 range), will be converted to original scale
- False: parameters are already in original scale, use as-is
Returns
-------
Union[np.array, tuple]
streamflow or (streamflow, states)
"""
model_param_dict = get_model_param_config("gr4j", kwargs)
param_ranges = model_param_dict["param_range"]
processed_params = process_parameters(
parameters, param_ranges, normalized=normalized_params
)
n_basins = processed_params.shape[0]
n_main = p_and_e.shape[0] - warmup_length
streamflow_ = np.empty((n_main, n_basins))
ets_out = np.empty((n_main, n_basins))
s_final = np.empty(n_basins)
r_final = np.empty(n_basins)
# Loop over basins in Python; each basin's hot path is the single fused
# JIT call below. For single-basin SCE-UA (the common case) this loop
# executes exactly once per evaluation.
for j in range(n_basins):
x1 = float(processed_params[j, 0])
x2 = float(processed_params[j, 1])
x3 = float(processed_params[j, 2])
x4 = float(processed_params[j, 3])
# Ensure contiguous float64 1D slices for numba (cheap if already so).
prcp = np.ascontiguousarray(p_and_e[:, j, 0], dtype=np.float64)
pet = np.ascontiguousarray(p_and_e[:, j, 1], dtype=np.float64)
if warmup_length > 0:
# Warm up by running the same core over the warmup window; we keep
# only the final storage states. No Python-level recursion into
# gr4j() — this used to double the work per evaluation.
_, _, s_warm, r_warm = _gr4j_jit_core(
prcp[:warmup_length],
pet[:warmup_length],
x1,
x2,
x3,
x4,
0.5 * x1,
0.5 * x3,
)
s_init, r_init = s_warm, r_warm
else:
s_init, r_init = 0.5 * x1, 0.5 * x3
qsim_j, ets_j, s_end, r_end = _gr4j_jit_core(
prcp[warmup_length:],
pet[warmup_length:],
x1,
x2,
x3,
x4,
s_init,
r_init,
)
streamflow_[:, j] = qsim_j
ets_out[:, j] = ets_j
s_final[j] = s_end
r_final[j] = r_end
streamflow = np.expand_dims(streamflow_, axis=2)
return (
(streamflow, ets_out, s_final, r_final)
if return_state
else (streamflow, ets_out)
)
production(p_and_e, x1, s_level=None)
¶
an one-step calculation for production store in GR4J the dimension of the cell: [batch, feature] Parameters
p_and_e P is pe[:, 0] and E is pe[:, 1]; similar with the "input" in the RNNCell
X1
Storage reservoir parameter;
s_level s_level means S in the GR4J Model; similar with the "hx" in the RNNCell Initial value of storage in the storage reservoir. Returns
tuple contains the Pr and updated S
Source code in hydromodel/models/gr4j.py
def production(
p_and_e: np.array, x1: np.array, s_level: Optional[np.array] = None
) -> Tuple[np.array, np.array]:
"""
an one-step calculation for production store in GR4J
the dimension of the cell: [batch, feature]
Parameters
----------
p_and_e
P is pe[:, 0] and E is pe[:, 1]; similar with the "input" in the RNNCell
x1:
Storage reservoir parameter;
s_level
s_level means S in the GR4J Model; similar with the "hx" in the RNNCell
Initial value of storage in the storage reservoir.
Returns
-------
tuple
contains the Pr and updated S
"""
# Calculate net precipitation and evapotranspiration
precip_difference = p_and_e[:, 0] - p_and_e[:, 1]
precip_net = np.maximum(precip_difference, 0.0)
evap_net = np.maximum(-precip_difference, 0.0)
if s_level is None:
s_level = 0.6 * x1
# s_level should not be larger than x1
s_level = np.clip(s_level, a_min=np.full(s_level.shape, 0.0), a_max=x1)
# Calculate the fraction of net precipitation that is stored
precip_store = calculate_precip_store(s_level, precip_net, x1)
# Calculate the amount of evaporation from storage
evap_store = calculate_evap_store(s_level, evap_net, x1)
# Update the storage by adding effective precipitation and
# removing evaporation
s_update = s_level - evap_store + precip_store
# s_level should not be larger than self.x1
s_update = np.clip(s_update, a_min=np.full(s_update.shape, 0.0), a_max=x1)
# Update the storage again to reflect percolation out of the store
perc = calculate_perc(s_update, x1)
s_update = s_update - perc
# perc is always lower than S because of the calculation itself, so we don't need clamp here anymore.
# The precip. for routing is the sum of the rainfall which
# did not make it to storage and the percolation from the store
current_runoff = perc + (precip_net - precip_store)
# TODO: check if evap_store is the real ET
return current_runoff, evap_store, s_update
routing(q9, q1, x2, x3, r_level=None)
¶
the GR4J routing-module unit cell for time-sequence loop Parameters
q9 q1 x2 Catchment water exchange parameter x3 Routing reservoir parameters r_level Beginning value of storage in the routing reservoir. Returns
Source code in hydromodel/models/gr4j.py
def routing(
q9: np.array, q1: np.array, x2, x3, r_level: Optional[np.array] = None
):
"""
the GR4J routing-module unit cell for time-sequence loop
Parameters
----------
q9
q1
x2
Catchment water exchange parameter
x3
Routing reservoir parameters
r_level
Beginning value of storage in the routing reservoir.
Returns
-------
"""
if r_level is None:
r_level = 0.7 * x3
# r_level should not be larger than self.x3
r_level = np.clip(r_level, a_min=np.full(r_level.shape, 0.0), a_max=x3)
groundwater_ex = x2 * (r_level / x3) ** 3.5
r_updated = np.maximum(
np.full(r_level.shape, 0.0), r_level + q9 + groundwater_ex
)
qr = r_updated * (1.0 - (1.0 + (r_updated / x3) ** 4) ** -0.25)
r_updated = r_updated - qr
qd = np.maximum(np.full(groundwater_ex.shape, 0.0), q1 + groundwater_ex)
q = qr + qd
return q, r_updated
s_curves1(t, x4)
¶
Unit hydrograph ordinates for UH1 derived from S-curves.
Source code in hydromodel/models/gr4j.py
@jit(nopython=True)
def s_curves1(t, x4):
"""
Unit hydrograph ordinates for UH1 derived from S-curves.
"""
if t <= 0:
return 0
elif t < x4:
return (t / x4) ** 2.5
else: # t >= x4
return 1
s_curves2(t, x4)
¶
Unit hydrograph ordinates for UH2 derived from S-curves.
Source code in hydromodel/models/gr4j.py
@jit(nopython=True)
def s_curves2(t, x4):
"""
Unit hydrograph ordinates for UH2 derived from S-curves.
"""
if t <= 0:
return 0
elif t < x4:
return 0.5 * (t / x4) ** 2.5
elif t < 2 * x4:
return 1 - 0.5 * (2 - t / x4) ** 2.5
else: # t >= x4
return 1
uh_gr4j(x4)
¶
Generate the convolution kernel for the convolution operation in routing module of GR4J
Parameters¶
x4 the dim of x4 is [batch] Returns
list UH1s and UH2s for all basins
Source code in hydromodel/models/gr4j.py
def uh_gr4j(x4):
"""
Generate the convolution kernel for the convolution operation in routing module of GR4J
Parameters
----------
x4
the dim of x4 is [batch]
Returns
-------
list
UH1s and UH2s for all basins
"""
uh1_ordinates = []
uh2_ordinates = []
for i in range(len(x4)):
n_uh1 = int(math.ceil(x4[i]))
n_uh2 = int(math.ceil(2.0 * x4[i]))
uh1_ordinate = np.zeros(n_uh1)
uh2_ordinate = np.zeros(n_uh2)
for t in range(1, n_uh1 + 1):
uh1_ordinate[t - 1] = s_curves1(t, x4[i]) - s_curves1(t - 1, x4[i])
for t in range(1, n_uh2 + 1):
uh2_ordinate[t - 1] = s_curves2(t, x4[i]) - s_curves2(t - 1, x4[i])
uh1_ordinates.append(uh1_ordinate)
uh2_ordinates.append(uh2_ordinate)
return uh1_ordinates, uh2_ordinates