This shows you the differences between two versions of the page.
| — |
research_pdmcode [2025/06/14 09:31] (current) trinh created |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| + | ====== An example of a conceptual rainfall-runoff model ====== | ||
| + | Here is an example implementation of the Probability Distributed Model. | ||
| + | |||
| + | <code matlab > | ||
| + | function [solution, hydrograph, obj] = simulate(obj, | ||
| + | | ||
| + | % Calculate length of each time step | ||
| + | dt = t_max / nt; | ||
| + | | ||
| + | % If precipitation is specified with a single value set precipitation | ||
| + | % rate to be the same for all time steps | ||
| + | if length(P) == 1 | ||
| + | P = P * ones(1, nt); | ||
| + | end | ||
| + | | ||
| + | % Initialise structure to save store values (an extra entry is | ||
| + | % included to store the initial values for each store) | ||
| + | | ||
| + | solution.C = zeros(1, nt+1); | ||
| + | solution.S = zeros(1, nt+1); | ||
| + | solution.Sf = zeros(1, nt+1); % fast storage value | ||
| + | solution.Ss = zeros(1, nt+1); % slow storage value | ||
| + | | ||
| + | % Update stores in every time step | ||
| + | for i = 1:nt+1 | ||
| + | | ||
| + | % Calculate the soil mosisture value based on the current soil | ||
| + | % moisture capacity | ||
| + | |||
| + | if obj.C < obj.par.c_max | ||
| + | S = obj.par.S_max * (1 - (1 - obj.C / obj.par.c_max) .^ (obj.par.b + 1)); | ||
| + | else | ||
| + | S = obj.par.S_max; | ||
| + | end | ||
| + | | ||
| + | % Save store values (for i=1 initial value is saved) | ||
| + | solution.C(i) = obj.C; | ||
| + | solution.S(i) = S; | ||
| + | solution.Sf(i) = obj.Sf; | ||
| + | solution.Ss(i) = obj.Ss; | ||
| + | | ||
| + | % If the last step was reached end the simulation | ||
| + | if i == nt + 1 | ||
| + | break | ||
| + | end | ||
| + | | ||
| + | % Calculate the drainage function | ||
| + | drainage = max(0, (S - obj.par.st) / obj.par.kg); | ||
| + | | ||
| + | % Calculate Pi (source term in ODE for the soil moisture storage) | ||
| + | Pi = P(i) - drainage; | ||
| + | | ||
| + | % Save old value of c (for runoff calculation) | ||
| + | C_old = obj.C; | ||
| + | | ||
| + | % Update value of soil moisture capacity, c | ||
| + | obj.C = min(obj.C + Pi * dt, obj.par.c_max); | ||
| + | | ||
| + | % Calculate the volume of the surface runoff from the water balance | ||
| + | runoff = Pi * dt - obj.par.S_max * ... | ||
| + | ((1 - C_old / obj.par.c_max) ^ (obj.par.b + 1) - ... | ||
| + | (1 - obj.C / obj.par.c_max) ^ (obj.par.b + 1)); | ||
| + | | ||
| + | % Calculare the rate of runoff generation | ||
| + | runoff = runoff / dt; | ||
| + | | ||
| + | % Calculate slow and fast flow | ||
| + | flow_s = exp(obj.Ss / obj.par.ks); | ||
| + | flow_f = obj.Sf / obj.par.kf; | ||
| + | | ||
| + | % Update values of the slow and fast storage | ||
| + | obj.Sf = obj.Sf + (runoff - flow_f) * dt; | ||
| + | obj.Ss = obj.Ss + (drainage - flow_s) * dt; | ||
| + | end | ||
| + | | ||
| + | % For each time step calculate the corresponding time (t), | ||
| + | % fast flow (Qf), slow flow (Qs) and total flow (Q_total) | ||
| + | hydrograph.t = linspace(0, t_max, nt+1); | ||
| + | hydrograph.Qf = solution.Sf / obj.par.kf; | ||
| + | hydrograph.Qs = exp(solution.Ss / obj.par.ks); | ||
| + | hydrograph.Q_total = hydrograph.Qs + hydrograph.Qf; | ||
| + | end | ||
| + | </ | ||