Trinh @ Bath

Trace:

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

— 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, 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
 +    </code>
Last modified: 2025/06/14 09:31