====== An example of a conceptual rainfall-runoff model ====== Here is an example implementation of the Probability Distributed Model. function [solution, hydrograph, obj] = simulate(obj, P, t_max, nt) % 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); % soil moisture storage capacity solution.S = zeros(1, nt+1); % soil moisture storage value 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