Initial commit

.gitattributes

@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto

.gitignore

@@ -0,0 +1,4 @@
+.*/
+*.pqo
+phreeqc.log
+*.fig

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2022 Matteo
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,56 @@
+# CWM1-matlab
+
+## A MATLAB implementation of Constructed Wetland Model 1 (CWM1) for reactive–transport simulators
+
+The Constructed Wetland Model 1 (CWM1) is a general model which describes the biochemical transformation/degradation processes for organic matter, nitrogen and sulphur in constructed wetlands (Langergraber et al. 2009 [[1]](#1)).
+
+This code is designed to be coupled with a chemical species transport simulator. Both for horizontal flow (HFCW) and vertical flow (VFCW) constructed wetlands the transport simulator usually consists of a variably saturated flow model (Richards' equation) and transport of particulates and solutes (advection-dispersion equation).
+
+The model includes the interpolation of some temperature-dependent parameters af a function of the water temperature with the method suggested by Henze et al. 2006 [[2]](#2).
+
+The CWM1 model consists of 17 processes (reactions) and 16 components (8 soluble and 8 particulate) which are represented by a system of ordinary differential equations (ODEs). The ODEs are solved using an explicit Runge-Kutta (2,3) method (`ode23` MATLAB function).
+
+## Usage
+
+CWM1-MATLAB basically consists of a main function `cwm1.m` which can be called as follows:
+
+```matlab
+C = cwm1(dt,params,init_cond)
+```
+
+The function takes the following inputs:
+
+- `dt`: simulation duration.
+- `params`: vector of model parameters. See `parameters.m` file.
+- `init_cond`: matrix [Nx16] of initial conditions, where N is the number of instances and 16 is the number of components of the CWM1 model
+
+As the function is meant to be coupled with a transport code (e.g., using a sequential non-iterative approach, SNIA), `dt` corresponds to the coupling step interval. `N` is the number of instances, i.e., the number of mesh nodes of the coupled transport model.
+
+The function returns a matrix of concentrations `C` of the final state of the `ode23` solver. Intermediate steps taken by the solver are discarded.
+
+## MEX function
+
+The `cwm1` is also made available as mex executables which are able to significantly increase the execution speed (> 5x faster than the .m file).
+Two versions are available in the release section:
+
+- `cwm1_mex.mexa64`: x86_64 Linux version, compiled with MATLAB R2019a on Ubuntu 18.04.4.
+- `cwm1_mex.mexw64`: x86_64 Windows version, compiled with MATLAB R2019b  on Windows 10.
+
+The compatibility with other systems is not guaranteed. The user is advised to build a specific executable with MATLAB Coder.
+
+## Benchmarks
+
+The file `main.m` shows an example of CWM1-MATLAB usage. The routine performs a speed comparison between plain MATLAB code and compiled .mex files.
+The routine also performs the validation of the code by comparison with CWM1 implemented in the PHREEQC software by Boog et al. 2018 [[3]](#3).
+
+![Benchmark](img/benchmark.png)
+
+The lines represent the simulated concentrations and circles represent the results from the PHREEQC model.
+
+## References
+
+<a id="1">[1]</a> Langergraber, G., Rousseau, D. P., García, J., & Mena, J. (2009). CWM1: a general model to describe biokinetic processes in subsurface flow constructed wetlands. *Water Science and Technology*, *59*(9), 1687-1697. DOI: [10.2166/wst.2009.131](https://doi.org/10.2166/wst.2009.131)
+
+<a id="2">[2]</a> Henze, M., Gujer, W., Mino, T., & van Loosdrecht, M. C. (2006). Activated Sludge Models ASM1, ASM2, ASM2d and ASM3. *IWA publishing*. ISBN (book): 9781900222242. DOI: [10.2166/9781780402369](https://doi.org/10.2166/9781780402369)
+
+<a id="3">[3]</a> Boog, J. (2018). Application: Treatment Wetlands. In *OpenGeoSys Tutorial* (pp. 63-90). SpringerBriefs in Earth System Sciences. Springer, Cham. DOI: [10.1007/978-3-319-67153-6_7](https://doi.org/10.1007/978-3-319-67153-6_7)

cwm1.m

@@ -0,0 +1,37 @@
+function C = cwm1(dt, params, init_cond)
+%% CWM1-MATLAB - HANDLING FUNCTION
+% CWM1 Initialize and run the ODE solver.
+%
+% Model parameters are passed as a vector to maximize compatibility with  
+% MATLAB Coder for mex building.
+%
+% Usage:
+%
+%   C = CWM1(dt,params,init_cond)
+%
+% Input:
+%   - dt: final time
+%   - params: vector of model parameters. See parameters.m file.
+%   - init_cond: matrix [Nx16], where N is the number of instances
+%     and 16 is the number of components of the CWM1 model
+%
+% Output:
+%   - C: Concentration matrix [Nx16] at the last integration step t = dt
+%
+% (c) Matteo M. 2022
+
+% Function handle system of equation
+ode_eqn = @(t,C) cwm1_odesystem(t,C,params); 
+
+C = zeros(size(init_cond)); % Initialize concentration matrix
+tspan = [0 dt];             % Time span
+
+for i = 1:size(init_cond,1)
+    % Solve (ode23 solver is faster than ode45, with acceptable accuracy)
+    [~,Ctmp] = ode23(ode_eqn, tspan, init_cond(i,:)); 
+
+    % Only save the final concentration at t = dt
+    C(i,:) = Ctmp(end,:); 
+end
+
+end

cwm1_odesystem.m

@@ -0,0 +1,149 @@
+function dCdt = cwm1_odesystem(t,C,p)
+%% CWM1-MATLAB - ODE SYSTEM
+% Definition of the system of ODEs
+%   N = 17 model processes
+%   M = 16 model components
+%
+% Equations from Langergraber et al. (2009)
+%
+% (c) Matteo M. 2022
+
+%% STOICHIOMETRIC MATRIX
+S = zeros(17,16);
+% comp#1 So
+S(2,1) = 1 - 1/p(49);
+S(4,1) = 1 - 1/p(49);
+S(7,1) = -(4.57-p(50))/p(50);
+S(15,1) = -(2-p(54))/p(54);
+% comp#2 Sf
+S(1,2) = 1 - p(46);
+S(2,2) = -1/p(49);
+S(3,2) = -1/p(49);
+S(6,2) = p(47);
+S(8,2) = p(47);
+S(9,2) = -1/p(51);
+S(10,2) = p(47);
+S(12,2) = p(47);
+S(14,2) = p(47);
+S(17,2) = p(47);
+% comp#3 Sa
+S(4,3) = -1/p(49);
+S(5,3) = -1/p(49);
+S(9,3) = (1-p(51))/p(51);
+S(11,3) = -1/p(52);
+S(13,3) = -1/p(53);
+% comp#4 Si
+S(1,4) = p(46);
+% comp#5 Snh
+S(1,5) = p(57) - (1-p(46))*p(55) - p(46)*p(56);
+S(2,5) = p(55)/p(49) - p(59);
+S(3,5) = S(2,5);
+S(4,5) = -p(59);
+S(5,5) = S(4,5);    
+S(6,5) = p(59) - p(47)*p(55) - (1 - p(47) - p(48))*p(57) - p(48)*p(58);
+S(7,5) = -p(59) - 1/p(50);
+S(8,5) = S(6,5);
+S(9,5) = p(55)/p(51) - p(59);
+S(10,5) = S(6,5);
+S(11,5) = S(4,5);
+S(12,5) = S(6,5);
+S(13,5) = S(4,5);
+S(14,5) = S(6,5);
+S(15,5) = S(4,5);
+S(16,5) = S(4,5);
+S(17,5) = S(6,5);
+% comp#6 Sno
+S(3,6) = -(1-p(49))/(2.86*p(49));
+S(5,6) = S(3,6);
+S(7,6) = 1/p(50);
+S(16,6) = -(1-p(54))/(0.875*p(54));
+% comp#7 Sso4
+S(13,7) = -(1-p(53))/(2*p(53));
+S(15,7) = 1/p(54);
+S(16,7) = 1/p(54);
+% comp#8 Sh2s
+S(13,8) = (1-p(53))/(2*p(53));
+S(15,8) = -1/p(54);
+S(16,8) = -1/p(54);
+% comp#9 Xs
+S(1,9) = -1;
+S(6,9) = 1 - p(47) - p(48);
+S(8,9) = S(6,9);
+S(10,9) = S(6,9);
+S(12,9) = S(6,9);
+S(14,9) = S(6,9);
+S(17,9) = S(6,9);
+% comp#10 Xi
+S(6,10) = p(48);
+S(8,10) = p(48);
+S(10,10) = p(48);
+S(12,10) = p(48);
+S(14,10) = p(48);
+S(17,10) = p(48);
+% comp#11 Xh
+S(2,11) = 1;
+S(3,11) = 1;
+S(4,11) = 1;
+S(5,11) = 1;
+S(6,11) = -1;
+% comp#12 Xa
+S(7,12) = 1;
+S(8,12) = -1;
+% comp#13 Xfb
+S(9,13) = 1;
+S(10,13) = -1;
+% comp#14 Xamb
+S(11,14) = 1;
+S(12,14) = -1;
+% comp#15 Xasrb
+S(13,15) = 1;
+S(14,15) = -1;
+% comp#16 Xsob
+S(15,16) = 1;
+S(16,16) = 1;
+S(17,16) = -1;
+
+
+%% PROCESS RATES
+P = zeros(1,17);
+% process#1 - Hydrolysis
+P(1) = p(1)* C(9) / (C(11)+C(13)) / ( p(2)+( C(9)/(C(11)+C(13)) )  ) * ( C(11)+p(3)*C(13)  );
+% process#2 - Aerobic growth of Xh on Sf
+P(2) = p(4) * C(2)/(p(8)+C(2)) * C(2)/(C(2)+C(3)) * C(1)/(p(7)+C(1)) * C(5)/(p(11)+C(5)) * p(12)/(p(12)+C(8)) * C(11);
+% process#3 - Anoxic growth of Xh on Sf
+P(3) = p(5)*p(4) * C(2)/(p(8)+C(2)) * C(2)/(C(2)+C(3)) * p(7)/(p(7)+C(1)) * C(6)/(p(10)+C(6)) * C(5)/(p(11)+C(5)) * p(12)/(p(12)+C(8)) * C(11);
+% process#4 - Aerobic growth of Xh on Sa
+P(4) = p(4) * C(3)/(p(9)+C(3)) * C(3)/(C(2)+C(3)) * C(1)/(p(7)+C(1)) * C(5)/(p(11)+C(5)) * p(12)/(p(12)+C(8)) * C(11);
+% process#5 - Anoxic growth of Xh on Sa
+P(5) = p(5)*p(4) * C(3)/(p(9)+C(3)) * C(3)/(C(2)+C(3)) * p(7)/(p(7)+C(1)) * C(6)/(p(10)+C(6)) * C(5)/(p(11)+C(5)) * p(12)/(p(12)+C(8)) * C(11);
+% process#6 - Lysis of Xh
+P(6) = p(6)*C(11);
+% process#7 - Aerobic growth of Xa on Snh
+P(7) = p(13) * C(5)/(p(16)+C(5)) * C(1)/(p(15)+C(1)) * p(17)/(p(17)+C(8)) * C(12);
+% process#8 - Lysis of Xa
+P(8) = p(14)*C(12);
+% process#9 - Growth of Xfb
+P(9) = p(18) * C(2)/(p(21)+C(2)) * p(24)/(p(24)+C(8)) * p(20)/(p(20)+C(1)) * p(22)/(p(22)+C(6)) * C(5)/(p(23)+C(5)) * C(13);
+% process#10 - Lysis of Xfb
+P(10) = p(19)*C(13);
+% process#11 - Growth of Xamb
+P(11) = p(25) * C(3)/(p(28)+C(3)) * p(31)/(p(31)+C(8)) * p(27)/(p(27)+C(1)) * p(29)/(p(29)+C(6)) * C(5)/(p(30)+C(5)) * C(14);
+% process#12 - Lysis of Xamb
+P(12) = p(26)*C(14);
+% process#13 - Growth of Xasrb
+P(13) = p(32) * C(3)/(p(35)+C(3)) * C(7)/(p(38)+C(7)) * p(39)/(p(39)+C(8)) * p(34)/(p(34)+C(1)) * p(36)/(p(36)+C(6)) * C(5)/(p(37)+C(5)) * C(15);
+% process#14 - Lysis of Xasrb
+P(14) = p(33)*C(15);
+% process#15 - Aerobic growth of Xsob on Sh2s
+P(15) = p(40) * C(8)/(p(45)+C(8)) * C(1)/(p(42)+C(1)) * C(5)/(p(44)+C(5)) * C(16);
+% process#16 - Anoxic growth of Xsob on Sh2s
+P(16) = p(40)* p(60) * C(8)/(p(45)+C(8)) * C(6)/(p(43)+C(6)) * p(42)/(p(42)+C(1)) * C(5)/(p(44)+C(5)) * C(16);
+% process#17 - Lysis of Xsob
+P(17) = p(41)*C(16);
+
+
+%% Build the system of ODEs
+dCdt = S'*P';
+
+
+end

main.m

@@ -0,0 +1,112 @@
+%% CWM1-MATLAB - Constructed Wetland Model 1
+% MATLAB implementation of CWM1 model.
+%
+% This code is designed to be coupled with a chemical species transport
+% simulator (e.g., water flow in unsaturated porous media through Richards'
+% Equation and an advection-dispersion model for solute transport)
+%
+% This file (main.m) shows an example of CWM1-MATLAB usage and performs a 
+% speed comparison between plain MATLAB code and compiled .mex files.
+% Moreover a benchmark is carried out to validate this model against the 
+% simulation results obtained with PHREEQC software. 
+%
+% References:
+%   CWM1 model: Langergraber et al. (2009)
+%   PHREEQC implementation: Boog (2018)
+%
+% (c) Matteo M. 2022
+clear;
+
+%% SIMULATION PARAMETERS
+N = 100;           % Number of instances (number of nodes of transport model)
+t_sim = 5;         % Simulation duration (days)
+dt = 0.02;         % Coupling time step (days)
+
+
+%% MODEL COMPONENTS
+% Initial conditions (g/L)
+% S = SOLUBLE
+So = 9.0e-3;            % comp#1. Dissolved oxygen
+Sf = 0.5e-3;            % comp#2. Fermentable, readily biodegradable COD
+Sa = 1.03e-3;			% comp#3. Fermentation products as acetate
+Sin	= 0.5e-3;			% comp#4. Inert soluble COD.
+Snh	= 0.5e-3;			% comp#5. Ammonium NH4+ and ammonia NH3 nitrogen
+Sno	= 40.0e-3;			% comp#6. Nitrate NO3- and nitrite NO2- nitrogen
+Sso4 = 10.00e-3;        % comp#7. Sulphate sulphur.
+Sh2s = 10.00e-3;		% comp#8. Dihydrogensulphide sulphur.
+% X = PARTICULATE
+Xs = 0.5e-3;			% comp#9. Slowly biodegradable particulate COD.
+Xi = 0.1e-3;			% comp#10. Inert particulate COD.
+Xh = 1.36e-4;           % comp#11. Heterotrophic bacteria.
+Xa = 1.36e-4;           % comp#12. Autotrophic nitrifying bacteria.
+Xfb = 1.36e-4;          % comp#13. Fermenting bacteria.
+Xamb = 1.36e-4;         % comp#14. Acetotrophic methanogenic bacteria.
+Xasrb = 1.36e-4;		% comp#15. Acetotrophic sulphate reducing bacteria.
+Xsob = 1.36e-4;         % comp#16. Sulphide oxidising bacteria.
+
+% Initial condition vector (mg/L)
+C0 = [So Sf Sa Sin Snh Sno Sso4 Sh2s Xs Xi Xh Xa Xfb Xamb Xasrb Xsob]*1000; 
+
+% Water temperature (°C)
+T = 20;
+
+%% MODEL PARAMETERS
+parameters;             % Load temperature-corrected model parameters from separate file
+
+M = length(C0);         % Number of component
+t = 0:dt:t_sim;         % Time vector
+nt = length(t);         % Number of time steps
+
+
+%% SOLVE - 1 - REGULAR MATLAB FUNCTION
+disp('_______Simulations started_______')
+disp(['Number of model instances: ' num2str(N)])
+init_cond = repmat(C0,N,1); 
+C = zeros(N,M,nt);      % Initialize state matrix
+C(:,:,1) = init_cond;   % Initial condition
+
+tic;
+for k = 2:nt
+    C(:,:,k) = cwm1(dt,params,C(:,:,k-1));
+end
+speed_reg = toc;
+disp(['Regular function, total elapsed time: ' num2str(speed_reg) ' s'])
+
+
+%% SOLVE - 2 - MEX FUCNTION SPEED TEST (>5x faster)
+tic;
+for k = 2:nt
+	C(:,:,k) = cwm1_mex(dt,params,C(:,:,k-1));
+end
+speed_mex = toc;
+disp(['MEX function, total elapsed time: ' num2str(speed_mex) ' s'])
+
+
+%% BENCKMARK COMPARISON VS PHREEQC RESULTS
+
+% LOAD PHREEQC OUTPUT
+filename = 'phreeqc/phout_sel.dat';
+opts = detectImportOptions(filename); % Preserves compatibility among different MATLAB versions
+phr_data = readtable(filename,opts);
+t_phr = phr_data.time/3600/24; t_phr(1) = 0;
+phr = table2array(phr_data(:,9:24));
+
+% PLOTS
+figure;
+plot(t,squeeze(C(1,:,:)),'Linewidth',1)
+lgnd = {'So','Sf','Sa','Sin','Snh','Sno','Sso4','Sh2s','Xs','Xi','Xh','Xa','Xfb','Xamb','Xasrb','Xsob'};
+hold all
+plot(t_phr,phr*1000, '.')
+legend(lgnd)
+xlabel('Elapsed time (days)')
+ylabel('Concentration (mg/L)')
+
+
+
+
+
+
+
+
+
+