Skip to content

GitLab

Commit c6bcd1b8 authored by Nicolas Rodriguez's avatar Nicolas Rodriguez
Browse files

Upload New File

parent 6f504fd2
Hide whitespace changes
Inline Side-by-side
model_calc.m 0 → 100644
View file @ c6bcd1b8
%
% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
%
% Author: Nicolas Rodriguez (nicolas.bjorn.rodriguez@gmail.com)
%
function [model_output,options]=model_calc(model_input,param,y0, options)
% solves the water balance, tracer balance, and master equations numerically
% Calculation period
calc=options.periods.calc;
lg=length(calc)-1; % size of the computation
% age options
dT=options.ages.dT; % age resolution for the calculations, in hours
save_interval=options.ages.save_interval; % interval of time (in the simulation) to save age calculations, in hours (usually a multiple of dT)
if length(calc)==length(options.periods.warmup)
if options.ages.save_warmup==0
save_interval=Inf;
end
end
init_RTD=options.ages.init_RTD; % initial residence time distribution in the storage unit
age_prctiles=options.ages.age_prctiles; % percentiles of the TTDs and RTDs to be calculated
T_group=options.ages.T_group; % groups to classify water younger than the given ages [days]
% Load data
t=model_input.t(calc);
dt=options.dt; % in hours
tau=dt/dT; % numerical dimensionless parameter, similar to a diffusion term. For example, tau=10% means that in dt hours, 10% of the water of a given age T has aged to T+dT
% tau is thus responsible for the diffusion-like profile of the TTDs.
tfull=model_input.tfull;
J=model_input.J(calc);
Cp=model_input.Cp(calc);
Cpfull=model_input.Cpfull;
PET=model_input.PET(calc);
Q=model_input.Q(calc);
size_true_warmup=model_input.size_true_warmup;
if length(calc)==length(options.periods.warmup)
ofst=max(0,size_true_warmup-length(options.periods.warmup)); % offset for convolution
else
ofst=max(size_true_warmup,length(options.periods.warmup)); % offset for convolution
end
% Get fixed parameters
for k=1:options.param.fix.nb
temp=options.param.fix.name{k};
eval(strcat(temp,'=',num2str(options.param.fix.val(k)),';'));
end
% Get calibration parameters
for k=1:options.param.cal.nb
temp=options.param.cal.name{k};
eval(strcat(temp,'=',num2str(param(k)),';'));
end
% Pre-allocation
S=zeros(lg+1,1);
M=zeros(lg+1,1);
C=zeros(lg+1,1);
C_S=zeros(lg+1,1);
C_Q=zeros(lg+1,1);
C_pest=zeros(lg+1,1);
C_weath=zeros(lg+1,1);
C_NO3=zeros(lg+1,1);
C_ET=zeros(lg+1,1);
ET=zeros(lg+1,1);
sum_ps=zeros(lg+1,1);
sum_pQ=zeros(lg+1,1);
sum_pET=zeros(lg+1,1);
% Initialization
if isempty(y0)
% Initial storages
S(1)=S_i;
% Initial isotopic masses
M(1)=S_i*C_i;
% Initial AVERAGE concentrations
C(1)=C_i;
else
S(1)=y0(1);
C(1)=y0(2);
M(1)=C(1)*S(1);
end
% Initial RTD
if isempty(init_RTD) % check if user specified initial RTDs
Ti=T0*365.25*24/dT; % estimated mean residence time of the catchment (in increments of dT)
szin=round(7*Ti); % size of the initialisation RTDs, here n=7 times the mean residence time (in increments of dT)
ps=1/(T0*365.25*24)*exp(-dT*(0:szin-1)/(T0*365.25*24));
ps=ps/(sum(ps)*dT);
else
ps=init_RTD;
end
% Pre-allocation
sum_ps(1)=sum(ps)*dT;
nb_columns=fix(lg/ceil(save_interval/dt)); % number of columns of the age storage matrices
store_index=zeros(1,nb_columns);
store_t=zeros(nb_columns,1);
store_ps=zeros(length(ps)+lg+1,nb_columns);
store_pQ=zeros(length(ps)+lg+1,nb_columns);
store_pET=zeros(length(ps)+lg+1,nb_columns);
store_wQ=zeros(length(ps)+lg,nb_columns);
store_wET=zeros(length(ps)+lg,nb_columns);
prct_S=zeros(nb_columns,length(age_prctiles)+1);
prct_Q=zeros(nb_columns,length(age_prctiles)+1);
prct_ET=zeros(nb_columns,length(age_prctiles)+1);
frc_S=zeros(length(T_group),nb_columns);
frc_Q=zeros(length(T_group),nb_columns);
frc_ET=zeros(length(T_group),nb_columns);
% pre-calculate the gamma distributions if needed
if options.ages.lookup==1
dS=options.ages.dS;
S_T=0:dS:2*S_i;
switch options.ages.SAStype_Q
case 'gam'
G_Q_1=gamcdf(S_T,SASQ_p1,SASQ_p2);
G_Q_2=[];
G_Q_3=[];
case '2g'
G_Q_1=gamcdf(S_T,SASQ_p3,SASQ_p5);
G_Q_2=gamcdf(S_T,SASQ_p4,SASQ_p6);
G_Q_3=[];
case '3g'
G_Q_1=gamcdf(S_T,SASQ_p4,SASQ_p7);
G_Q_2=gamcdf(S_T,SASQ_p5,SASQ_p8);
G_Q_3=gamcdf(S_T,SASQ_p6,SASQ_p9);
otherwise
G_Q_1=[]; G_Q_2=[]; G_Q_3=[];
end
switch options.ages.SAStype_ET
case 'gam'
G_ET_1=gamcdf(S_T,SASET_p1,SASET_p2);
G_ET_2=[];
G_ET_3=[];
case '2g'
G_ET_1=gamcdf(S_T,SASET_p3,SASET_p5);
G_ET_2=gamcdf(S_T,SASET_p4,SASET_p6);
G_ET_3=[];
case '3g'
G_ET_1=gamcdf(S_T,SASET_p4,SASET_p7);
G_ET_2=gamcdf(S_T,SASET_p5,SASET_p8);
G_ET_3=gamcdf(S_T,SASET_p6,SASET_p9);
otherwise
G_ET_1=[]; G_ET_2=[]; G_ET_3=[];
end
else
S_T=[];
G_Q_1=[]; G_Q_2=[]; G_Q_3=[];
G_ET_1=[]; G_ET_2=[]; G_ET_3=[];
end
% solver options
if strcmp(options.solver.type,'IE')
tol=options.solver.tol;
prec=options.solver.prec;
maxiter=options.solver.maxiter;
dis=options.solver.display;
max_disp=options.solver.maxdisp;
ndisp=0;
end
% Data saving options
duration=0;
save_index=1;
nchecks=1;
MTTD_Q=0;
MTTD_ET=0;
MRTD=0;
msg_S=0;
% Main loop
for j=1:lg
% Check if not too many messages have been displayed
if strcmp(options.solver.type,'IE')
if strcmp(dis,'on')
if ndisp>max_disp
dis='off';
warning('solver: STOPPING ERROR MESSAGE DISPLAY BECAUSE TOO MANY RECURRING MESSAGES!');
end
end
end
% Regular checks on the solutions
if j==ceil(options.solver.checks*nchecks*lg)
nchecks=nchecks+1;
% Checks on state variables
y=cat(2,S,M,C);
cond=['p';'n';'n']; % conditions on the signs of elements of y: pos (storage), neg (isotope mass), neg (isotope concentration)
[feas,errtypes]=checksol(y,cond);
if feas==0
model_output.feas=feas;
model_output.errtypes=errtypes;
return;
end
end
% Water Balance calculations
if strcmp(options.solver.type,'IE')
% implicit Euler
eq=@(x)(x-S(j)-J(j+1)+dt*(Q(j)+Psi*PET(j)*tanh((x/S_ET)^m)));
% AET decreases from PET towards 0 when S<S_ET, using a smoothing parameter m
lb=S(j)-100;
ub=S(j)+100;
% numerical solution to the equation by method of bisection (slowest but safest)
[S(j+1),~,~,~,exitflag]=mysolver(eq,lb,ub,tol,prec,maxiter,dis);
if exitflag<2
ndisp=ndisp+1;
end
ET(j)=Psi*PET(j)*tanh((S(j+1)/S_ET)^m);
else % explicit Euler
ET(j)=Psi*PET(j)*tanh((S(j)/S_ET)^m);
% Global smoother equivalent to max(S_iter,0), without the numerically-induced instabilities
S_iter=S(j)+J(j+1)-dt*(Q(j)+ET(j));
S(j+1)=0.5*(S_iter+sqrt(S_iter^2+n));
if msg_S==0
dS_rel=abs(S(j+1)-S_iter)/S_iter*100;
if dS_rel>10 % 10% difference with the smoother?
msg_S=1;
end
end
end
% Tracer Mass Balance calculations
% calculate SAS functions and TTDs
[p_Q,omega_Q]=evalSAS(dT,ps,S(j),0,options.ages.SAStype_Q,SASQ_p1,SASQ_p2,SASQ_p3,SASQ_p4,SASQ_p5,SASQ_p6,SASQ_p7,SASQ_p8,SASQ_p9,S_T,G_Q_1,G_Q_2,G_Q_3);
[p_ET,omega_ET]=evalSAS(dT,ps,S(j),0,options.ages.SAStype_ET,SASET_p1,SASET_p2,SASET_p3,SASET_p4,SASET_p5,SASET_p6,SASET_p7,SASET_p8,SASET_p9,S_T,G_ET_1,G_ET_2,G_ET_3);
% pesticide input
Cp_pest=zeros(size(Cpfull)); % 0 mg/L pesticide all the time
Cp_pest(end-round(8*365.25*24/dt))=50; % 50 mg/L during the 100 mm precip event on '01-Oct-2009 00:00'
Cp_N=60*ones(size(Cpfull)); % 60 mg/L, constant nitrate input
% decrease of the nitrate input from 60 mg/L to 0 mg/L in 150 days, starting on '01-Oct-2009 00:00'
Cp_N(end-round(8*365.25*24/dt):end-round(8*365.25*24/dt)+150)=60:-60/150:0;
Cp_N(end-round(8*365.25*24/dt)+151:end)=0;
% interpolate data vectors (necessary for the convolution operation)
t_int=(tfull(j+ofst):(-dT/24):tfull(1))';
if length(t_int)==1
Cp_int=Cpfull(j+ofst);
Cp_pest_int=Cp_pest(j+ofst);
Cp_N_int=Cp_N(j+ofst);
else
temp=interp1(tfull(1:j+ofst),Cpfull(1:j+ofst),flip(t_int),'nearest');
Cp_int=reshape(temp,[length(temp) 1]);
temp=interp1(tfull(1:j+ofst),Cp_pest(1:j+ofst),flip(t_int),'nearest');
Cp_pest_int=reshape(temp,[length(temp) 1]);
temp=interp1(tfull(1:j+ofst),Cp_N(1:j+ofst),flip(t_int),'nearest');
Cp_N_int=reshape(temp,[length(temp) 1]);
end
% prepare the vectors for convolution
ps_conv=ps;
if length(ps)>=length(Cp_int)
Cp_conv=cat(1,flip(Cp_int),zeros(length(ps)-length(Cp_int),1));
Cp_pest_conv=cat(1,flip(Cp_pest_int),zeros(length(ps)-length(Cp_pest_int),1));
Cp_N_conv=cat(1,flip(Cp_N_int),zeros(length(ps)-length(Cp_N_int),1));
else
Cp_conv=flip(Cp_int);
Cp_pest_conv=flip(Cp_pest_int);
Cp_N_conv=flip(Cp_N_int);
ps_conv(end:length(Cp_int))=0;
p_Q(end:length(Cp_int))=0;
p_ET(end:length(Cp_int))=0;
end
% account for pesticide decay
decay_factor=1/(24*50); % 1/hr, here it corresponds to tau_dec=50 days
[r,~]=size(Cp_pest_conv);
% account for weathering products
Ceq=200; % equilibirum concentration in mg/L
k=1/(24*550); % solute kinetic constant 1/hr, here it corresponds to tau_eq=550 days
C_prod=Ceq*(1-exp(-k*dT*(0:1:r-1)))'; % kinetic law for increasing solute concentration towards Ceq
Cp_pest_conv=Cp_pest_conv.*transpose(exp(-decay_factor*dT*(0:1:r-1)));
% do the convolution
C_S(j)=dT*ps_conv*Cp_conv;
C_Q(j)=dT*p_Q*Cp_conv;
C_pest(j)=dT*p_Q*Cp_pest_conv;
C_weath(j)=dT*p_Q*C_prod;
C_NO3(j)=dT*p_Q*Cp_N_conv;
C_ET(j)=dT*p_ET*Cp_conv;
% update the tracer mass balance
M(j+1)=M(j)+J(j+1)*Cp(j+1)-dt*(Q(j)*C_Q(j)+ET(j)*C_ET(j));
C(j+1)=M(j+1)/S(j+1);
% Age calculations
% Replace NaN with 0 for calculations
om_Q=omega_Q;
om_Q(isnan(om_Q))=0;
om_ET=omega_ET;
om_ET(isnan(om_ET))=0;
% numerical scheme:
% forward 1st order finite difference in time
% backward 1st order finite difference in age
% recombined in a vectorized version
a=(S(j)*(1-tau)-dt*(Q(j)*om_Q+ET(j)*om_ET))/S(j+1);
b=tau*S(j)/S(j+1);
c=J(j+1)/S(j+1);
% [ REMOVAL ] [ AGING ] [ NEW WATER ]
ps_iter=cat(2,a,0).*cat(2,ps,0)+b*cat(2,0,ps(1:end))+c*cat(2,1/dT,0*ps);
% Neglect small numbers in the vectors
ps_iter(ps_iter<eps)=0;
% shorten vectors
if options.ages.short==1
z=find(ps_iter==0);
if ~isempty(z)
% make sure the vectors are not shortened too much
cp=cumsum(ps_iter*dT);
u=find(cp<0.95);
if ~isempty(u)
i=min(z(z>=u(end)));
if ~isempty(i)
ps_iter=ps_iter(1:i);
end
end
end
end
% update RTDs
ps=ps_iter;
% check if distributions sum to 1
sum_ps(j+1)=sum(ps)*dT;
sum_pQ(j)=sum(p_Q)*dT;
sum_pET(j)=sum(p_ET)*dT;
% normalize distributions to 1 if asked
if options.ages.norm==1
ps=ps/sum_ps(j+1);
end
% Calculate Flow weighted Master Travel Time / Residence Time Distributions
tailsize=length(p_Q)-length(MTTD_Q); % difference of size between previous mttd and the TTD to be added
if tailsize>0
MTTD_Q(end:length(p_Q))=0; % resizing of the previous mttd to avoid errors when adding the current TTD
MTTD_Q=MTTD_Q+Q(j)*dt*p_Q;
else
MTTD_Q(1:length(p_Q))=MTTD_Q(1:length(p_Q))+Q(j)*dt*p_Q; % else no need to resize
end
tailsize=length(p_ET)-length(MTTD_ET); % difference of size between previous mttd and the TTD to be added
if tailsize>0
MTTD_ET(end:length(p_ET))=0; % resizing of the previous mttd to avoid errors when adding the current TTD
MTTD_ET=MTTD_ET+ET(j)*dt*p_ET;
else
MTTD_ET(1:length(p_ET))=MTTD_ET(1:length(p_ET))+ET(j)*dt*p_ET; % else no need to resize
end
tailsize=length(ps)-length(MRTD);
if tailsize>0
MRTD(end:length(ps))=0;
MRTD=MRTD+S(j)*ps;
else
MRTD(1:length(ps))=MRTD(1:length(ps))+S(j)*ps;
end
% save results
duration=duration+dt;
if ~isinf(save_interval)
if duration>=save_interval
duration=0; % reset
store_index(save_index)=j;
store_t(save_index)=t(j);
store_ps(1:length(ps),save_index)=ps/(sum(ps)*dT);
store_pQ(1:length(p_Q),save_index)=p_Q/(sum(p_Q)*dT);
store_pET(1:length(p_ET),save_index)=p_ET/(sum(p_ET)*dT);
store_wQ(1:length(omega_Q),save_index)=omega_Q;
store_wET(1:length(omega_ET),save_index)=omega_ET;
% Calculate age metrics
[prct_S(save_index,:),frc_S(:,save_index)]=age_metrics(ps,dT,T_group,age_prctiles);
[prct_Q(save_index,:),frc_Q(:,save_index)]=age_metrics(p_Q,dT,T_group,age_prctiles);
[prct_ET(save_index,:),frc_ET(:,save_index)]=age_metrics(p_ET,dT,T_group,age_prctiles);
save_index=save_index+1; % selection of the next column of the storage matrices
end
end
end
% last time step: (same steps as during the main loop)
ET(j+1)=Psi*PET(j)*tanh((S(j+1)/S_ET)^m);
[p_Q,~]=evalSAS(dT,ps,S(j+1),0,options.ages.SAStype_Q,SASQ_p1,SASQ_p2,SASQ_p3,SASQ_p4,SASQ_p5,SASQ_p6,SASQ_p7,SASQ_p8,SASQ_p9,S_T,G_Q_1,G_Q_2,G_Q_3);
[p_ET,~]=evalSAS(dT,ps,S(j+1),0,options.ages.SAStype_ET,SASET_p1,SASET_p2,SASET_p3,SASET_p4,SASET_p5,SASET_p6,SASET_p7,SASET_p8,SASET_p9,S_T,G_ET_1,G_ET_2,G_ET_3);
sum_pQ(j+1)=sum(p_Q)*dT;
sum_pET(j+1)=sum(p_ET)*dT;
t_int=(tfull(j+ofst):(-dT/24):tfull(1))';
if length(t_int)==1
Cp_int=Cpfull(j+ofst+1);
else
temp=interp1(flip(tfull(j+ofst+1:-1:1)),flip(Cpfull(j+ofst+1:-1:1)),flip(t_int),'nearest');
Cp_int=reshape(temp,[length(temp) 1]);
end
ps_conv=ps;
if length(ps)>=length(Cp_int)
Cp_conv=cat(1,flip(Cp_int),zeros(length(ps)-length(Cp_int),1));
else
Cp_conv=flip(Cp_int);
ps_conv(end:length(Cp_int))=0;
p_Q(end:length(Cp_int))=0;
p_ET(end:length(Cp_int))=0;
end
C_S(j+1)=dT*ps_conv*Cp_conv;
C_Q(j+1)=dT*p_Q*Cp_conv;
C_ET(j+1)=dT*p_ET*Cp_conv;
% Checks state variables for NaN, Inf...
y=cat(2,S,M,C);
cond=['p';'n';'n']; % conditions on the signs of elements of y
[feas,errtypes]=checksol(y,cond);
% Check numerical precision
num_wat_bal=(S(end)-(S(1)+sum(J(2:end))-sum((Q(1:end-1)+ET(1:end-1))*dt)))/mean(S)*100; % numerical water balance
num_tra_bal=(S(end)*C(end)-(S(1)*C(1)+sum(J(2:end).*Cp(2:end))-sum((Q(1:end-1).*C_Q(1:end-1)+ET(1:end-1).*C_ET(1:end-1))*dt)))/mean(M)*100; % numerical tracer balance
if abs(num_wat_bal)>5 || abs(num_tra_bal)>5
feas=0;
errtypes=cat(2,errtypes,5);
end
% Check numerical age precision
if mean(sum_ps)<0.9 || mean(sum_ps)>1.1 || prctile(sum_ps,25)<0.9 || prctile(sum_ps,75)>1.1
feas=0;
errtypes=cat(2,errtypes,6);
end
if options.display==1
if msg_S==1
warning('model_calc: Q had to be unaccounted for at some point to prevent negative S values');
end
end
% replace 0 in C by NaN
C(C==0)=NaN;
C_S(C_S==0)=NaN;
C_Q(C_Q==0)=NaN;
C_ET(C_ET==0)=NaN;
% normalisation of the MTTDs by total outflow volume
MTTD_Q=MTTD_Q./sum(Q(1:lg)*dt);
MTTD_Q=1/sum(MTTD_Q*dT)*MTTD_Q;
MTTD_ET=MTTD_ET./sum(ET(1:lg)*dt);
MTTD_ET=1/sum(MTTD_ET*dT)*MTTD_ET;
MRTD=MRTD./sum(S(1:lg));
MRTD=1/sum(MRTD*dT)*MRTD;
% Flux weighted age quantities
[FYW_Q,MTT_Q,mTT_Q]=weighted_age_metrics(MTTD_Q,frc_Q,dT,Q(store_index));
[FYW_ET,MTT_ET,mTT_ET]=weighted_age_metrics(MTTD_ET,frc_ET,dT,ET(store_index));
[FYW_S,MRT,mRT]=weighted_age_metrics(MRTD,frc_S,dT,S(store_index));
% Smooth the simulated stream/ET concentration
if options.movavg==1
C_S=movmean(C_S,round(dT/dt));
C_Q=movmean(C_Q,round(dT/dt));
C_ET=movmean(C_ET,round(dT/dt));
C_pest=movmean(C_pest,round(dT/dt));
C_weath=movmean(C_weath,round(dT/dt));
C_NO3=movmean(C_NO3,round(dT/dt));
end
% export results
model_output.ages.sum_ps=sum_ps;
model_output.ages.sum_pQ=sum_pQ;
model_output.ages.sum_pET=sum_pET;
if ~isempty(store_ps)
model_output.ages.store_ps=store_ps;
else
model_output.ages.store_ps=ps';
end
model_output.ages.store_pQ=store_pQ;
model_output.ages.store_pET=store_pET;
model_output.ages.store_wQ=store_wQ;
model_output.ages.store_wET=store_wET;
model_output.ages.store_t=store_t;
model_output.ages.MRTD=MRTD;
model_output.ages.MTTD_Q=MTTD_Q;
model_output.ages.MTTD_ET=MTTD_ET;
model_output.ages.prct_S=prct_S;
model_output.ages.prct_Q=prct_Q;
model_output.ages.prct_ET=prct_ET;
model_output.ages.frc_S=frc_S;
model_output.ages.frc_Q=frc_Q;
model_output.ages.frc_ET=frc_ET;
model_output.ages.FYW_Q=FYW_Q;
model_output.ages.FYW_ET=FYW_ET;
model_output.ages.FYW_S=FYW_S;
model_output.ages.MTT_Q=MTT_Q;
model_output.ages.MTT_ET=MTT_ET;
model_output.ages.MRT=MRT;
model_output.ages.mTT_Q=mTT_Q;
model_output.ages.mTT_ET=mTT_ET;
model_output.ages.mRT=mRT;
model_output.C(calc)=C;
model_output.C_S(calc)=C_S;
model_output.S(calc)=S;
model_output.M(calc)=M;
model_output.C_Q(calc)=C_Q;
model_output.C_pest(calc)=C_pest;
model_output.C_NO3(calc)=C_NO3;
model_output.C_weath(calc)=C_weath;
model_output.C_ET(calc)=C_ET;
model_output.ET(calc)=ET;
model_output.num_wat_bal=num_wat_bal;
model_output.num_tra_bal=num_tra_bal;
model_output.feas=feas;
model_output.errtypes=errtypes;
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Luxembourg Institute of Science and Technology - www.list.lu.