MARRMOT_23
Implementation
MARRMoT v2.1.1 rev eeb7e15 m_23_lascam_24p_3s; standardized continuous structure, not original model
Time step: daily
Backend: octave-cli + Octave optim
Calibrated parameters: 24
Temperature required: no
Parameters and initial configuration
Parameter |
Supported calibration range |
Default |
|---|---|---|
|
0 to 200 |
100.0 |
|
0 to 5 |
2.5 |
|
1 to 2000 |
1000.5 |
|
0.01 to 0.99 |
0.5 |
|
0.01 to 0.99 |
0.5 |
|
0.01 to 0.99 |
0.5 |
|
0 to 5 |
2.5 |
|
0 to 10 |
5.0 |
|
0 to 5 |
2.5 |
|
0 to 10 |
5.0 |
|
0 to 200 |
100.0 |
|
0 to 5 |
2.5 |
|
0 to 1 |
0.5 |
|
0 to 1 |
0.5 |
|
0 to 10 |
5.0 |
|
0 to 1 |
0.5 |
|
0 to 1 |
0.5 |
|
0.01 to 200 |
100.005 |
|
0 to 1 |
0.5 |
|
0 to 10 |
5.0 |
|
0.01 to 200 |
100.005 |
|
1 to 5 |
3.0 |
|
0 to 1 |
0.5 |
|
0 to 10 |
5.0 |
|
Fixed initial/configuration value |
0.0 |
|
Fixed initial/configuration value |
0.0 |
|
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;
af = theta(1); % Catchment-scale infiltration parameter [mm/d]
bf = theta(2); % Catchment-scale infiltration non-linearity parameter [-]
ac = theta(7); % Variable contributing area scaling [-]
bc = theta(8); % Variable contributing area non-linearity [-]
ass = theta(9); % Subsurface saturation area scaling [-]
bss = theta(10); % Subsurface saturation area non-linearity [-]
c = theta(11); % Maximum infiltration rate [mm/d]
ag = theta(12); % Interception base parameter [mm/d]
bg = theta(13); % Interception fraction parameter [-]
gf = theta(14); % F-store evaporation scaling [-]
df = theta(15); % F-store evaporation non-linearity [-]
td = theta(16); % Recharge time parameter [d-1]
ab = theta(17); % Groundwater flow scaling [-]
bb = theta(18); % Groundwater flow base rate [mm/d]
ga = theta(19); % A-store evaporation scaling [-]
da = theta(20); % A-store evaporation non-linearity [-]
aa = theta(21); % Subsurface storm flow rate [mm/d]
ba = theta(22); % Subsurface storm flow non-linearity [-]
gb = theta(23); % B-store evaporation scaling [-]
db = theta(24); % B-store evaporation non-linearity [-]
% auxiliary parameters
aux_theta = obj.aux_theta;
amax = aux_theta(1); % Maximum contributing area depth [mm]
fmax = aux_theta(2); % Infiltration depth scaling [mm]
bmax = aux_theta(3); % Groundwater depth scaling [mm]
amin = aux_theta(4); % Minimum contributing area depth [mm]
% delta_t
delta_t = obj.delta_t;
% stores
S1 = S(1);
S2 = S(2);
S3 = S(3);
% 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
tmp_phiss= area_1(ass,bss,S2,amin,amax);
tmp_phic = area_1(ac,bc,S2,amin,amax);
tmp_fss = infiltration_5(af,bf,S3,bmax,S1,fmax);
flux_pg = interception_5(bg,ag,P);
flux_ei = effective_1(P,flux_pg);
flux_qse = saturation_11(ac,bc,S2,amin,amax,flux_pg);
flux_pc = infiltration_4(flux_pg-flux_qse,c);
flux_qie = effective_1(flux_pg-flux_qse,flux_pc);
flux_qsse= saturation_12(tmp_phiss,tmp_phic,flux_pc);
flux_fa = infiltration_4(max(0,flux_pc*min(1,(1-tmp_phiss)/(1-tmp_phic))),tmp_fss);
flux_qsie= effective_1(flux_pc,flux_fa+flux_qsse);
flux_ef = evap_19(gf,df,S1,fmax,Ep,delta_t);
flux_rf = recharge_3(td,S1);
flux_ea1 = evap_1(S2,tmp_phic*Ep,delta_t) ;
flux_ea2 = evap_19(ga,da,S2,amax,Ep,delta_t);
flux_qa = saturation_11(aa,ba,S2,amin,amax,1);
flux_ra = recharge_4(tmp_phic,tmp_fss,delta_t);
flux_qb = baseflow_8(bb,ab,S3,bmax);
flux_eb = evap_19(gb,db,S3,bmax,Ep,delta_t);
% stores ODEs
dS1 = flux_fa - flux_ef - flux_rf;
dS2 = flux_qsse + flux_qsie + flux_qb - flux_ea1 - ...
flux_ea2 - flux_ra - flux_qa;
dS3 = flux_rf + flux_ra - flux_eb - flux_qb;
% outputs
dS = [dS1 dS2 dS3];
fluxes = [flux_ei, flux_pg, flux_qse, flux_qie,...
flux_pc, flux_qsse, flux_qsie, flux_fa,...
flux_ef, flux_rf, flux_ea1, flux_ea2,...
flux_qa, flux_ra, flux_qb, flux_eb];
end
% STEP runs at the end of every timestep
function [out] = area_1(p1,p2,S,Smin,Smax,varargin)
%area_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: Auxiliary function that calculates a variable contributing area.
% Constraints: A <= 1
% @(Inputs): p1 - linear scaling parameter [-]
% p2 - exponential scaling parameter [-]
% S - current storage [mm]
% Smin - minimum contributing storage [mm]
% Smax - maximum contributing storage [mm]
% varargin(1) - smoothing variable r (default 0.01)
% varargin(2) - smoothing variable e (default 5.00)
if size(varargin,2) == 0
out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin)); % default smoothing
elseif size(varargin,2) == 1
out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin,varargin(1))); % user-specified smoothing
elseif size(varargin,2) == 2
out = min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin,varargin(1),varargin(2))); % user-specified smoothing
end
end
function [out] = infiltration_5(p1,p2,S1,S1max,S2,S2max)
%infiltration_5
% 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: Maximum infiltration rate non-linearly based on relative deficit and storage
% Constraints: S2 >= 0 - prevents complex numbers
% f <= 10^9 - prevents numerical issues with Inf outcomes
% @(Inputs): p1 - base infiltration rate [mm/d]
% p2 - exponential scaling parameter [-]
% S1 - current storage in S1 [mm]
% S1max - maximum storage in S1 [mm]
% S2 - current storage in S2 [mm]
% S2max - maximum storage in S2 [mm]
out = max(0,min(10^9,p1.*(1-S1./S1max).*max(0,S2./S2max).^(-1.*p2)));
end
function [out] = interception_5(p1,p2,In)
%interception_5
% 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 after a combined absolute amount and fraction are intercepted
% Constraints: f >= 0
% @(Inputs): p1 - fraction that is not throughfall [-]
% p2 - constnat interception and evaporation [mm/d]
% In - incoming flux [mm/d]
out = max(p1.*In-p2,0);
end
function [out] = effective_1(In1,In2)
%effective_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: General effective flow (returns flux [mm/d])
% Constraints: In1 > In2
% @(Inputs): In1 - first flux [mm/d]
% In2 - second flux [mm/d]
out = max(In1-In2,0);
end
function [out] = saturation_11(p1,p2,S,Smin,Smax,In,varargin)
%saturation_11
% 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 flow from a store with different degrees
% of saturation (min exponential variant)
% Constraints: f <= In
% @(Inputs): p1 - linear scaling parameter [-]
% p2 - exponential scaling parameter [-]
% S - current storage [mm]
% Smin - minimum contributing storage [mm]
% Smax - maximum contributing storage [mm]
% In - incoming flux [mm/d]
% varargin(1) - smoothing variable r (default 0.01)
% varargin(2) - smoothing variable e (default 5.00)
if size(varargin,2) == 0
out = In.*min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin));
elseif size(varargin,2) == 1
out = In.*min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin,varargin(1)));
elseif size(varargin,2) == 2
out = In.*min(1,p1.*(max(0,S-Smin)./(Smax-Smin)).^p2).*...
(1-smoothThreshold_storage_logistic(S,Smin,varargin(1),varargin(2)));
end
end
function [out] = infiltration_4(fluxIn,p1)
%infiltration_4
% 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: Constant infiltration rate
% Constraints: f <= fin
% @(Inputs): p1 - Infiltration rate [mm/d]
% fin - incoming flux [mm/d]
out = min(fluxIn,p1);
end
function [out] = saturation_12(p1,p2,In)
%saturation_12
% 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 flow from a store with different degrees
% of saturation (min-max linear variant)
% Constraints: f >= 0
% @(Inputs): p1 - maximum contributing fraction area [-]
% p2 - minimum contributing fraction area [-]
% In - incoming flux [mm/d]
out = max(0,(p1-p2)./(1-p2)).*In;
end
function [out] = evap_19(p1,p2,S,Smax,Ep,dt)
%evap_19 Creates scaled, non-linear evaporation
% 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: Non-linear scaled evaporation
% Constraints: f <= Ep
% f <= S/dt
% @(Inputs): p1 - linear scaling parameter [-]
% p2 - exponential scaling parameter [-]
% S - current storage [mm]
% Smax - maximum storage [mm]
% Ep - potential evapotranspiration rate [mm/d]
% dt - time step size [d]
out = min([S/dt,Ep,p1.*max(0,S/Smax).^(p2).*Ep]);
end
function [out] = recharge_3(p1,S)
%recharge_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: Linear recharge
% Constraints: -
% @(Inputs): p1 - time coefficient [d-1]
% S - current storage [mm]
out = p1.*S;
end
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] = recharge_4(p1,S,dt)
%recharge_4
% 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: Constant recharge
% Constraints: f <= S/dt
% @(Inputs): p1 - time coefficient [d-1]
% S - current storage [mm]
% dt - time step size [d]
out = min(p1,S/dt);
end
function [out] = baseflow_8(p1,p2,S,Smax)
%baseflow_8
% 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 scaled outflow from a deficit store
% Constraints: S <= Smax
% @(Inputs): p1 - base outflow rate [mm/d]
% p2 - exponential scaling parameter [-]
% S - current storage [mm]
% Smax - maximum contributing storage [mm]
out = p1.*(exp(p2.*min(1,max(S,0)./Smax))-1);
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_23").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.