MARRMOT_36

Implementation

MARRMoT v2.1.1 rev eeb7e15 m_36_modhydrolog_15p_5s; standardized continuous structure, not original model

  • Time step: daily

  • Backend: octave-cli + Octave optim

  • Calibrated parameters: 15

  • Temperature required: no

Inspect the exact implementation.

Parameters and initial configuration

Parameter

Supported calibration range

Default

p01

0 to 5

2.5

p02

0 to 600

300.0

p03

0 to 15

7.5

p04

1 to 2000

1000.5

p05

0 to 1

0.5

p06

0 to 1

0.5

p07

0 to 20

10.0

p08

0 to 50

25.0

p09

0 to 1

0.5

p10

0.99 to 1

0.995

p11

0 to 0.5

0.25

p12

-10 to 10

0.0

p13

0 to 1

0.5

p14

0 to 1

0.5

p15

0 to 100

50.0

s01

Fixed initial/configuration value

0.0

s02

Fixed initial/configuration value

0.0

s03

Fixed initial/configuration value

0.0

s04

Fixed initial/configuration value

0.0

s05

Fixed initial/configuration value

0.0

Ranges/defaults are implementation contracts, not universal priors or a recommended basin calibration. Consult source comments for parameter units and coupling.

Governing equations

The following source is the exact model kernel used by this adapter. For MARRMoT it includes the state derivative and each referenced flux function; the solver and routing are described above. Original notices and source citations are retained in the files.

        function [dS, fluxes] = model_fun(obj, S)
            % parameters
            theta = obj.theta;
            insc    = theta(1);     % Maximum interception capacity, [mm]
            coeff   = theta(2);     % Maximum infiltration loss parameter, [-]
            sq      = theta(3);     % Infiltration loss exponent, [-]
            smsc    = theta(4);     % Maximum soil moisture capacity, [mm]
            sub     = theta(5);     % Proportionality constant, [-]
            crak    = theta(6);     % Proportionality constant, [-]
            em      = theta(7);     % Plant-controled maximum evaporation rate [mm/d]
            dsc     = theta(8);     % Maximum depression capacity, [mm]
            ads     = theta(9);     % Land fraction functioning as depression storage, [-]
            md      = theta(10);    % Depression storage parameter, [-], default = 1
            vcond   = theta(11);    % Leakage coefficient, [mm/d]
            dlev    = theta(12);    % Datum around which groundwater fluctuates relative to river bed, [mm]
            k1      = theta(13);    % Flow exchange parameter, [d-1] 
            k2      = theta(14);    % Flow exchange parameter, [d-1] 
            k3      = theta(15);    % Flow exchange parameter, [d-1] 
            
            % delta_t
            delta_t = obj.delta_t;
            
            % stores
            S1 = S(1);
            S2 = S(2);
            S3 = S(3);
            S4 = S(4);
            S5 = S(5);
            
            % climate input
            t = obj.t;                             % this time step
            climate_in = obj.input_climate(t,:);   % climate at this step
            P  = climate_in(1);
            Ep = climate_in(2);
            T  = climate_in(3);
            
            % fluxes functions
            flux_Ei     = evap_1(S1,Ep,delta_t);         
            flux_EXC    = interception_1(P,S1,insc);
            flux_INF    = infiltration_1(coeff,sq,S2,smsc,flux_EXC);
            flux_INT    = interflow_1(sub,S2,smsc,flux_INF);
            flux_REC    = recharge_1(crak,S2,smsc,flux_INF-flux_INT);
            flux_SMF    = flux_INF - flux_INT - flux_REC;
            flux_Et     = evap_2(em,S2,smsc,Ep,delta_t);
            flux_GWF    = saturation_1(flux_SMF,S2,smsc);
            flux_RUN    = flux_EXC - flux_INF;
            flux_TRAP   = depression_1(ads,md,S3,dsc,flux_RUN,delta_t);
            flux_Ed     = evap_1(S3,ads*Ep,delta_t);
            flux_DINF   = ads .* infiltration_2(coeff,sq,S2,smsc,flux_SMF,S3,delta_t);
            flux_SEEP   = exchange_3(vcond,S4,dlev);
            flux_SRUN   = flux_RUN - flux_TRAP;
            flux_FLOW   = exchange_1(k1,k2,k3,S4,flux_SRUN,delta_t);
            flux_Q      = baseflow_1(1,S5);

            % stores ODEs
            dS1 = P          - flux_Ei    - flux_EXC;
            dS2 = flux_SMF   + flux_DINF  - flux_Et    - flux_GWF;
            dS3 = flux_TRAP  - flux_Ed    - flux_DINF;
            dS4 = flux_REC   + flux_GWF   - flux_SEEP  - flux_FLOW;
            dS5 = flux_SRUN  + flux_INT   + flux_FLOW  - flux_Q; 
            
            % outputs
            dS = [dS1 dS2 dS3 dS4 dS5];
            fluxes = [flux_Ei,   flux_EXC, flux_INF,  flux_INT,...
                      flux_REC,  flux_SMF, flux_Et,   flux_GWF,...
                      flux_TRAP, flux_Ed,  flux_DINF, flux_SEEP,...
                      flux_FLOW, flux_Q,   flux_RUN,  flux_SRUN];
        end
        
        % STEP runs at the end of every timestep
function [out] = evap_1(S,Ep,dt)
%evap_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Evaporation at the potential rate
% Constraints:  f <= S/dt
% @(Inputs):    S    - current storage [mm]
%               Ep   - potential evaporation rate [mm/d]
%               dt   - time step size

out = min(S/dt,Ep);

end
function [out] = interception_1(In,S,Smax,varargin)
%interception_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Interception excess when maximum capacity is reached
% Constraints:  -
% @(Inputs):    In   - incoming flux [mm/d]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               varargin(1) - smoothing variable r (default 0.01)
%               varargin(2) - smoothing variable e (default 5.00)

if size(varargin,2) == 0
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax));
elseif size(varargin,2) == 1
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax,varargin(1)));
elseif size(varargin,2) == 2
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax,varargin(1),varargin(2)));    
end

end
function [out] = infiltration_1(p1,p2,S,Smax,fin)
%infiltration_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Infiltration as exponentially declining based on relative storage
% Constraints:  f <= fin
% @(Inputs):    p1   - maximum infiltration rate [mm,/d]
%               p2   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               fin  - size of incoming flux [mm/d]

out = min(p1.*exp((-1*p2*S)./Smax),fin);

end
function [out] = interflow_1(p1,S,Smax,flux)
%interflow_1 interflow based on incoming flux size

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Interflow as a scaled fraction of an incoming flux
% Constraints:  -
% @(Inputs):    p1   - linear scaling parameter [-]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               flux - incoming flux [mm/d]

out = p1*S/Smax*flux;

end
function [out] = recharge_1(p1,S,Smax,flux)
%recharge_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Recharge as scaled fraction of incoming flux
% Constraints:  -
% @(Inputs):    p1   - fraction of flux that is recharge [-]
%               S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               flux - incoming flux [mm/d]

out = p1*S/Smax*flux;

end
function [out] = evap_2(p1,S,Smax,Ep,dt)
%evap_2 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Evaporation at a scaled, plant-controlled rate
% Constraints:  f <= Ep
%               f <= S/dt
% @(Inputs):    p1   - plant-controlled base evaporation rate [mm/d]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               Ep   - potential evapotranspiration rate [mm/d]
%               dt   - time step size [d]

out = min([p1*S/Smax,Ep,S/dt]);

end
function [out] = saturation_1(In,S,Smax,varargin)
%saturation_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Saturation excess from a store that has reached maximum capacity
% Constraints:  -
% @(Inputs):    In   - incoming flux [mm/d]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               varargin(1) - smoothing variable r (default 0.01)
%               varargin(2) - smoothing variable e (default 5.00)

if size(varargin,2) == 0
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax));
elseif size(varargin,2) == 1
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax,varargin(1)));
elseif size(varargin,2) == 2
    out = In.*(1-smoothThreshold_storage_logistic(S,Smax,varargin(1),varargin(2)));    
end

end
function [out] = depression_1(p1,p2,S,Smax,flux,dt)
%depression_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Exponential inflow to surface depression store
% Constraints:  f <= (Smax-S)/dt
%               S <= Smax
% @(Inputs):    p1   - linear scaling parameter [-]
%               p2   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               dt   - time step size [d]

out = min(p1.*exp(-1.*p2.*S./max(Smax-S,0)).*flux,max((Smax-S)/dt,0));

end
function [out] = infiltration_2(p1,p2,S1,S1max,flux,S2,dt)
%infiltration_2 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Infiltration as exponentially declining based on relative storage
% Constraints:  0 <= f <= S2/dt
% @(Inputs):    p1    - maximum infiltration rate [mm,/d]
%               p2    - exponential scaling parameter [-]
%               S1    - current storage in S1 [mm]
%               S1max - maximum storage in S1 [mm]
%               flux  - reduction of infiltration rate by infiltration 
%                       demand already fulfilled elsewhere [mm/d]
%               S2    - storage available for infiltration [mm]
%               dt    - time step size [d]

out = max(min((p1.*exp(-1*p2*S1./S1max))-flux,S2/dt),0);

end
function [out] = exchange_3(p1,S,p2)
%exchange_3 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Water exchange with infinite size store based on threshold
% Constraints:  -
% @(Inputs):    p1   - base leakage time delay [d-1]
%               p2   - threshold for flow reversal [mm]
%               S    - current storage [mm]

out = p1*(S-p2);

end
function [out] = exchange_1(p1,p2,p3,S,fmax,dt)
%exchange_1 two-way channel exchange: linear and exponential.

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Water exchange between aquifer and channel
% Constraints:  f <= fIn
% @(Inputs):    p1   - linear scaling parameter [-]
%               p2   - linear scaling parameter [-]
%               p3   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               fmax - maximum flux size [mm/d]
%               dt   - time step size [d]

out = max((p1*abs(S/dt) + p2*(1-exp(-1*p3*abs(S/dt)))).*sign(S),-1*fmax);

end
function [out] = baseflow_1(p1,S)
% baseflow_1 

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Outflow from a linear reservoir
% Constraints:  -
% @(Inputs):    p1   - time scale parameter [d-1]
%               S    - current storage [mm]

out = p1.*S;

end
function [out] = smoothThreshold_storage_logistic(S,Smax,r,e)
%smoothThreshold_storage_logistic Logisitic smoother for storage threshold functions.

% Copyright (C) 2018 Wouter J.M. Knoben
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

%   Smooths the transition of threshold functions of the form:
%
%   Q = { P, if S = Smax
%       { 0, if S < Smax
%
%   By transforming the equation above to Q = f(P,S,Smax,e,r):
%   Q = P * 1/ (1+exp((S-Smax+r*e*Smax)/(r*Smax)))
%
%   Inputs:
%   S       : current storage
%   Smax    : maximum storage
%   r       : [optional] smoothing parameter rho, default = 0.01
%   e       : [optional] smoothing parameter e, default 5
%
%   NOTE: this function only outputs the multiplier. This needs to be
%   applied to the proper flux utside of this function.
%
%   NOTE: can be applied for temperature thresholds as well (i.e. snow
%   modules). This simply means that S becomes T, and Smax T0.

% Check for inputs and use defaults if not provided
% NOTE: this is not very elegant, but it is more than a factor 10 faster then: 
% if ~exist('r','var'); r = 0.01; end
% if ~exist('e','var'); e = 5.00; end
if nargin == 2
    r = 0.01;
    e = 5.00;
elseif nargin == 3
    r = r{1};
    e = 5.00;
elseif nargin == 4
    r = r{1};
    e = e{1};
end

% Calculate multiplier
Smax = max(Smax,0);   % this avoids numerical instabilities when Smax<0
if r*Smax == 0
    out = 1 ./ (1+exp((S-Smax+r*e*Smax)/(r)));
else
    out = 1 ./ (1+exp((S-Smax+r*e*Smax)/(r*Smax)));
end

end

Note

This standardized MARRMoT structure is not identical to the original named model. p01... follow exact source order; s01... are fixed initial stores, defaulting to zero. Octave + optim are required. Solver behavior can differ across runtime versions.

Simulation

from basinforge import Basin, get_model

basin = Basin.from_csv("basin.csv", basin_id="A", area_km2=1200,
                       q_unit="m3/s", timestep="daily")
q_mm = get_model("MARRMOT_36").simulate(basin)
q_m3s = basin.to_m3s(q_mm)

Supply your actual data and catchment area; temperature-dependent models require a temperature column. For non-daily models, choose an appropriate warmup in model steps.

See calibration, input requirements, sources and verification limitations.