Description of the GR4J model

The GR4J model has four parameters to optimise during calibration:

X1: the production store maximal capacity (mm),
X2: the catchment water exchange coefficient (mm/day),
X3: the one-day maximal capacity of the routing reservoir (mm),
X4: the HU1 unit hydrograph time base (days).

We denote by P (mm/day) the rainfall amount and by E (mm/day) the potential evapotranspiration (PET). P is an estimation of the catchment rainfall and E can come from a mean interannual PET curve.
The following equations correspond to the equations integrated over a time step.
The first operation is the neutralization of P by E to determine a net rainfall Pn and a net evapotranspiration En, calculated by:
If P > E, then Pn = P – E and En = 0
If P < E, then Pn = 0 and En = E – P
In case Pn is different from zero, a fraction Ps of Pn goes to the production reservoir and is calculated by:

where X1 (mm) and S are, respectively, the maximum capacity and the production store level.
Otherwise, when En is different from zero, a part of evaporation Es is removed from the production store. It is given by:

The production store level is updated through:
S = S – Es + Ps
A percolation called Perc coming from the production store is then calculated:

The production store level is then again updated:
S = S – Perc
The water quantity Pr that finally reaches the routing part of the model is:
Pr = Perc + (Pn – Ps)
Pr is divided into two flow components, 90 % being routed by a unit hydrograph HU1 and a routing store and 10 % by a unique unit hydrograph HU2.
HU1 and HU2 depend on the same parameter X4, the time base of HU1 expressed in days.
The hydrograms ordinates are calculated from the S curves (the accumulation of the proportion of unit rainfall treated by the hydrogram in function of time), respectively named SH1 and SH2.
SH1 is defined in function of time by:
For t = 0
For 0 < t < X4
For t > X4
SH2 is defined in function of time by:
For t = 0
For 0 < t < X4
For X4 < t < 2X4
For t > 2X4
The ordinates of HU1 and HU2 are then obtained from:

where j is an integer.
For each time step i, the outputs Q9 and Q1 of the two hydrograms are calculated with:

with l = int(X4)+1 and m = int(2.X4)+1, with int(.) the integer part.
A groundwater exchange term (loss or gain) is calculated with:

with R the routing store level, X3 the one-day maximal capacity of the store and X2 the water exchange coefficient, which is positive in case of a gain, and negative in case of a loss, or null.
The level in the routing store is updated by adding the Q9 output of the hydrogram HU1 and F:
R = max (0 ; R + Q9 + F)
Then, it empties in an output Qr given by:

The level in the store becomes: R = R – Qr
The output Q1 of the hydrogram HU2 goes through the same exchanges to give the flow component Qd:
Qd = max (0 ; Q1+F)
The total streamflow Q is finally given by: Q = Qr + Qd
To know more: check our publications.