MARRMOT_33
Implementation
MARRMoT v2.1.1 rev eeb7e15 m_33_sacramento_11p_5s; standardized continuous structure, not original model
Time step: daily
Backend: octave-cli + Octave optim
Calibrated parameters: 11
Temperature required: no
Parameters and initial configuration
Parameter |
Supported calibration range |
Default |
|---|---|---|
|
0 to 1 |
0.5 |
|
1 to 2000 |
1000.5 |
|
0.005 to 0.995 |
0.5 |
|
0.005 to 0.995 |
0.5 |
|
0 to 1 |
0.5 |
|
0 to 7 |
3.5 |
|
0.005 to 0.995 |
0.5 |
|
0.005 to 0.995 |
0.5 |
|
0 to 1 |
0.5 |
|
0 to 1 |
0.5 |
|
0 to 1 |
0.5 |
|
Fixed initial/configuration value |
0.0 |
|
Fixed initial/configuration value |
0.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;
pctim = theta(1); % Fraction impervious area [-]
kuz = theta(5); % Interflow runoff coefficient [d-1]
rexp = theta(6); % Base percolation rate non-linearity factor [-]
pfree = theta(9); % Fraction of percolation directed to free water stores [-]
klzp = theta(10); % Primary baseflow runoff coefficient [d-1]
klzs = theta(11); % Supplemental baseflow runoff coefficient [d-1]
theta_d = obj.theta_derived;
uztwm = theta_d(1); % Maximum upper zone tension water storage [mm]
uzfwm = theta_d(2); % Maximum upper zone free water storage [mm]
lztwm = theta_d(3); % Maximum lower zone tension water storage [mm]
lzfwpm = theta_d(4); % Maximum lower zone primary free water storage [mm]
lzfwsm = theta_d(5); % Maximum lower zone supplemental free water storage [mm]
pbase = theta_d(6); % Base percolation rate [mm/d]
zperc = theta_d(7); % Base percolation rate multiplication factor [-]
% 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);
% fluxes functions
% Original formulation using MARRMoT fluxes is very slow on sacramento,
% individual functions have been explicitly coded underneath.
flux_qdir = split_1(pctim,P);
flux_peff = split_1(1-pctim,P);
flux_ru = soilmoisture_1(S1,uztwm,S2,uzfwm);
flux_euztw = evap_7(S1,uztwm,Ep,delta_t);
flux_twexu = saturation_1(flux_peff,S1,uztwm);
flux_qsur = saturation_1(flux_twexu,S2,uzfwm);
flux_qint = interflow_5(kuz,S2);
flux_euzfw = evap_1(S2,max(0,Ep-flux_euztw),delta_t);
flux_pc = percolation_4(pbase,zperc,rexp,max(0,lztwm-S3)+max(0,lzfwpm-S4)+max(0,lzfwsm-S5),lztwm+lzfwpm+lzfwsm,S2,uzfwm,delta_t);
flux_pctw = split_1(1-pfree,flux_pc);
flux_elztw = evap_7(S3,lztwm,max(0,Ep-flux_euztw-flux_euzfw),delta_t);
flux_twexl = saturation_1(flux_pctw,S3,lztwm);
flux_twexlp = split_1(deficitBasedDistribution(S4,lzfwpm,S5,lzfwsm),flux_twexl);
flux_twexls = split_1(deficitBasedDistribution(S5,lzfwsm,S4,lzfwpm),flux_twexl);
flux_pcfwp = split_1(pfree*deficitBasedDistribution(S4,lzfwpm,S5,lzfwsm),flux_pc);
flux_pcfws = split_1(pfree*deficitBasedDistribution(S5,lzfwsm,S4,lzfwpm),flux_pc);
flux_rlp = soilmoisture_2(S3,lztwm,S4,lzfwpm,S5,lzfwsm);
flux_rls = soilmoisture_2(S3,lztwm,S5,lzfwsm,S4,lzfwpm);
flux_qbfp = baseflow_1(klzp,S4);
flux_qbfs = baseflow_1(klzs,S5);
% stores ODEs
dS1 = flux_peff + flux_ru - flux_euztw - flux_twexu;
dS2 = flux_twexu - flux_euzfw - flux_qsur - flux_qint - flux_ru - flux_pc;
dS3 = flux_pctw + flux_rlp + flux_rls - flux_elztw - flux_twexl;
dS4 = flux_twexlp + flux_pcfwp - flux_rlp - flux_qbfp;
dS5 = flux_twexls + flux_pcfws - flux_rls - flux_qbfs;
% outputs
dS = [dS1, dS2, dS3, dS4, dS5];
fluxes = [flux_qdir, flux_peff, flux_ru, flux_euztw, flux_twexu,...
flux_qsur, flux_qint, flux_euzfw, flux_pc, flux_pctw,...
flux_elztw, flux_twexl, flux_twexlp, flux_twexls, flux_pcfwp,...
flux_pcfws, flux_rlp, flux_rls, flux_qbfp, flux_qbfs];
end
% STEP runs at the end of every timestep.
function [out] = split_1(p1,In)
%split_1 flow splitting
% 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: Split flow (returns flux [mm/d])
% Constraints: -
% @(Inputs): p1 - fraction of flux to be diverted [-]
% In - incoming flux [mm/d]
out = p1.*In;
end
function [out] = soilmoisture_1(S1,S1max,S2,S2max)
%soilmoisture_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: Water rebalance to equal relative storage (2 stores)
% Constraints: -
% @(Inputs): 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 = ((S2.*S1max-S1.*S2max)/(S1max+S2max)).*smoothThreshold_storage_logistic(S1./S1max,S2./S2max);
end
function [out] = evap_7(S,Smax,Ep,dt)
%evap_7 evaporation based on scaled current water storage, limited by
%potential rate.
% 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 scaled by relative storage
% Constraints: f <= S/dt
% @(Inputs): S - current storage [mm]
% Smax - maximum contributing storage [mm]
% Ep - potential evapotranspiration rate [mm/d]
% dt - time step size [d]
out = min(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] = interflow_5(p1,S)
%interflow_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: Linear interflow
% 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] = percolation_4(p1,p2,p3,p4,p5,S,Smax,dt)
%percolation_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: Demand-based percolation scaled by available moisture
% Constraints: f <= S/dt
% f >= 0 prevents erratic numerical behaviour
% @(Inputs): p1 - base percolation rate [mm/d]
% p2 - percolation rate increase due moisture deficiencies [mm/d]
% p3 - non-linearity parameter [-]
% p4 - summed deficiency across all model stores [mm]
% p5 - summed capacity of model stores [mm]
% S - current storage in the supplying store [mm]
% Smax - maximum storage in the supplying store [mm]
% dt - time step size [d]
% Note: for certain extreme parameter values (very small stores, highly
% non-linear p3) and small computational errors that lead to small negative
% S values, this function behaves erratically. The max(0,S/Smax) part
% prevents this. Similarly, the first max(0,...) part prevents negative
% percolation demands as a result of small numerical errors.
out = max(0,min(S/dt,max(0,S./Smax).*(p1.*(1+p2.*(p4./p5).^(1+p3)))));
end
function [ f1,f2 ] = deficitBasedDistribution( S1,S1max,S2,S2max )
%DEFICITBASEDDISTRIBUTION Calculates a fractional split for two stores,
%based on the relative deficit in each. Currently used in:
% m_33_sacramento_11p_5s
% 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.
%%CALCULATE RELATIVE DEFICITS
rd1 = (S1-S1max)/S1max;
rd2 = (S2-S2max)/S2max;
%%CALCULATE FRACTIONAL SPLIT
if rd1+rd2 ~= 0
% Deficit exists and can be used to compute the split
f1 = rd1/(rd1+rd2);
f2 = rd2/(rd1+rd2);
else
% Both deficits are zero, and we revert to distribution based on
% relative maximum store size
f1 = S1max/(S1max+S2max);
f2 = S2max/(S1max+S2max);
end
end
function [out] = soilmoisture_2(S1,S1max,S2,S2max,S3,S3max)
%soilmoisture_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: Water rebalance to equal relative storage (3 stores)
% Constraints: -
% @(Inputs): S1 - current storage in S1 [mm]
% S1max - maximum storage in S1 [mm]
% S2 - current storage in S2 [mm]
% S2max - maximum storage in S2 [mm]
% S3 - current storage in S3 [mm]
% S3max - maximum storage in S3 [mm]
out = (S2.*(S1.*(S2max+S3max)+S1max.*(S2+S3))./((S2max+S3max).*(S1max+S2max+S3max))).* ...
smoothThreshold_storage_logistic(S1./S1max,(S2+S3)./(S2max+S3max));
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_33").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.