In the environmental risk assessment of pesticides, required for the authorization of plant protection products, information is needed on the fate of such a pesticide in surface water. To that end, the long-term concentration of a substance in surface water, including organic matter in the water column, sediment layers, and flow, after single or multiple exposure of the surface water is simulated. One model that is extensively used for these simulations in Europe is TOXSWA. TOXSWA takes as user-definable substance-specific input the degradation rates of the substance in the water and sediment compartments (expressed as first order half lives DT50), and the organic matter-water partition coefficient (Kom). Usually, these substance-specific parameters are based on the results of laboratory experiments, more specifically the soil adsorption/desorption study (Kom) and the water/sediment degradation study (DT50s). One of the drawbacks of this approach is that in water/sediment studies the concentration of the substance over time in the water and sediment compartment is recorded. The concentration-vs-time curves for both compartments are therefore compound curves, reflecting both the degradation of the substance and the distribution over the separate compartments of the system. The concentration of the substance in the aqueous compartment e.g. is determined by the degradation of the substance in the water layer as well as the rate of adsorption to the sediment and possibly also the test vessel wall, and sometimes also the rate of volatilization. Figure 1 shows a schematic overview of the different processes that play a role in the fate of a compound in a water/sediment system.
This makes it fairly difficult to determine the 'real' degradation rates of a compound. Water/sediment study results usually report the overall 'dissipation' rate constants (or DT50s) from the different compartments. The recommended approach for performing TOXSWA simulations based on these dissipation rates therefore is to use the overall ('system') dissipation rate as a rough surrogate for the degradation rate in water, and to set the degradation DT50 in sediment to 10000 days. This provides a worst case estimate of the behaviour of the substance in water. However, statistical methods exist that can estimate the individual rate constants for systems as depicted above, if sufficient concentration-vs-time information for the separate compartments is available. According to Schrap et al. (1994), the fate of a substance in a water/sediment system as formalized in Figure 1 can be described by the following set of linked differential equations:
with Cw = concentration of the substance in water (mol L-1), Cs = concentration of the substance in sediment (mol kg-1), Cg = 'concentration' of the substance on the glass vessel wall (mol kg-1), S = ratio of sediment (as organic matter) to water (kg L-1), and G = ratio of glass to water (kg L-1). kx are first order rate constants (d-1), and KL is a Langmuir constant (L kg-1). As and added benefit, the organic matter/water partition coefficient KOM can be derived from the ratio of the system independent rate constants for exchange between water and organic matter (see also Schrap et al. (1994)):
Sets of linked differential equations can be solved numerically, even if no analytical solutions exist, using so-called ODE solvers for Initial Value Problems. If sufficient experimental data for the appearance and disappearance of the individual species in the different compartments of such a system are available, it is possible to estimate the individual rate constants by fitting the numerically solved set of differential equations to the experimental data, using a non-linear least squares algorithm. This can be done in a variety of mathematical modelling languages, or programs, such as ACSL, MATLAB, ModelMaker, Berkeley Madonna, to name a few. Due to historical reasons, I use MATLAB to achieve this. MATLAB is a general purpose numerical mathematics system, that offers, among other things, statistical, optimization and ODE routines. I will present here a sample set of MATLAB program files that can be used to perform the parameter estimation described above. In this example, it is assumed that all metabolic transformation steps can adequately be represented by first order non-equilibrium reactions, and that all intercompartment movement can be represented by first order equilibrium processes.
function dydt = water_der(t,y,pars) %substance in water; degradation to metabolite_1; %metabolite in water; degradation; %linear sorption of substance from water to sediment; %substance in sediment; degradation to metabolite_2; %metabolite_2 in sediment; degradation k_c_w = pars(1); k_m1_w = pars(2); k_c_wts = pars(3); k_c_stw = pars(4); k_c_s = pars(5); k_m2_s = pars(6); f_m1 = pars(7); f_m2 = pars(8); dydt = zeros(size(y)); C_c_w = y(1); C_m1_w = y(2); C_c_s_pseudo = y(3); C_m2_s_pseudo = y(4); dydt(1) = -(k_c_w + k_c_stw) * C_c_w + k_c_wts * C_c_s_pseudo; dydt(2) = f_m1 * k_c_w * C_c_w - k_m1_w * C_m1_w; dydt(3) = k_c_wts * C_c_w - (k_c_s + k_c_stw) * C_c_s_pseudo; dydt(4) = f_m2 * k_c_s * C_c_s_pseudo - k_m2_s * C_m2_s_pseudo;
function [tijd, C_c_w, C_m1_w, C_c_s_pseudo, C_m2_s_pseudo] = water(pars,stop_time,Y0) period = [0 stop_time]; tijd = []; C_c_w = []; C_m1_w = []; C_c_s_pseudo = []; C_m2_s_pseudo = []; [t,y] = ode15s(@water_der,period,Y0,[],pars); tijd = [tijd;t]; C_c_w = [C_c_w;y(:,1)]; C_m1_w = [C_m1_w;y(:,2)]; C_c_s_pseudo = [C_c_s_pseudo;y(:,3)]; C_m2_s_pseudo = [C_m2_s_pseudo;y(:,4)]; laatst = length (tijd);
function diff = fit_water(pars,X,Y,span,Y0) [tijd, Q_c_w,Q_m1_w,Q_c_s,Q_m2_s] = water(pars, span, Y0); diff = interp1(tijd, [Q_c_w,Q_m1_w,Q_c_s,Q_m2_s],X)-Y;
%Fitting of rate constants for a system with two compartments (water and
%sediment), possibly an adsorptive vessel wall, and degradation in the water and sediment
%compartments.
%For substance, degradation in the aqueous phase is to metabolite_1,
%in the sediment phase to metabolite_2
%See also S.M. Schrap et al., ENVIRONMENTAL SCIENCE AND POLLUTION
%RESEARCH INTERNATIONAL; 1 (1). 1994, pp. 21-28.
%Discrete time points
%X = [0;1;3;7;10;14;21;28;42];
%span = X(length(X));
if ~exist('SYSTEM')
SYSTEM = 0;
end
switch SYSTEM
case (1)
%Observed amounts of substance in water (Q_c_w), and sediment (Q_c_s), metabolite_1 in water (Q_m1_w) and metabolite_2 in sediment (Q_m2_s)
TITLE = 'Example';
X = [0;1;3;7;10;14;21;28;42];
span = X(length(X));
Y1 = [66.2;51.8;42.3;19.7;7.6;3.4;1.9;.2;0];
Y2 = [4.4;4.7;3.1;9;15.8;13.8;9.8;7.8;7.1];
Y3 = [21.1;22.3;17.8;10.1;5.9;4.8;3.1;.9;3];
Y4 = [1.5;3.4;5;9.2;10.5;6.7;5.9;3;1];
case (2)
otherwise
disp('Houston, we have a problem.');
return;
end
%initial values
Y0 = [70;0;0;0];
switch SYSTEM
case (1)
%Amount of water and sediment
A_ww = .075; %amount of water in L
A_ss = .0061; %amount of sediment in kg
F_om = 0.038; %fraction of organic matter in sediment;
otherwise
disp('Houston, we have a problem.');
return;
end
A_omss = F_om*A_ss; %amount of organic matter in system
om_to_water_ratio = A_omss/A_ww; %self-evident...
%Concentrations of compounds in water and sediment
C_c_w = Y1./A_ww;
C_m1_w = Y2./A_ww;
C_c_s = Y3./A_omss;
C_m2_s = Y4./A_omss;
Y0 = Y0./[A_ww;A_ww;A_omss;A_omss];
%pseudoconcentration of compounds in sediment
C_c_s_pseudo = C_c_s .* om_to_water_ratio;
C_m2_s_pseudo = C_m2_s .* om_to_water_ratio;
Y0 = Y0.*[1;1;om_to_water_ratio;om_to_water_ratio];
%initial guesses for rate constants and formation efficiencies
K0 = [.02 .01 .0001 .1 .01 .0001 .1 .1];
%Lower and Upper bounds for the rate constants
LB = [0 0 0 0 0 0 0 0];
UB = [Inf Inf Inf Inf Inf Inf 1 1];
Y = [C_c_w C_m1_w C_c_s_pseudo C_m2_s_pseudo];
options = [];
[est,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN] = lsqnonlin(@fit_water,K0,LB,UB,options,X,Y,span,Y0);
[tijd, C_c_w_est, C_m1_w_est, C_c_s_pseudo_est, C_m2_s_pseudo_est] = water(est, span, Y0);
C_c_s_est = C_c_s_pseudo_est ./ om_to_water_ratio;
C_m2_s_est = C_m2_s_pseudo_est ./ om_to_water_ratio;
plot(tijd,C_c_w_est,'b',tijd,C_m1_w_est,'g',tijd,C_c_s_pseudo_est,'r',tijd,C_m2_s_pseudo_est,'k',X,Y(:,1),'+b',X,Y(:,2),'og',X,Y(:,3),'dr',X,Y(:,4),'xk');
cf = nlparci (est,RESIDUAL,JACOBIAN);
kaas = [est' cf];
logKom = log10((kaas(3,1)/om_to_water_ratio)/kaas(4,1));
DT50_c_in_water = log(2)/kaas(1,1);
DT50_m1_in_water = log(2)/kaas(2,1);
DT50_c_in_sediment = log(2)/kaas(5,1);
DT50_m2_in_sediment = log(2)/kaas(6,1);
f_sub_to_m1_w = kaas(7,1);
f_sub_to_m2_s = kaas(8,1);
Half_max_water = interp1(C_c_w_est,tijd,C_c_w_est(1,1)/2);
tbl = [{'Param' 'Value' 'Lower bound' 'Upper bound'}];
tbl = [tbl;{'DT50 substance w (d)' DT50_c_in_water log(2)/kaas(1,3) log(2)/kaas(1,2)}];
tbl = [tbl;{'DT50 metabolite_1 w (d)' DT50_m1_in_water log(2)/kaas(2,3) log(2)/kaas(2,2)}];
tbl = [tbl;{'DT50 substance s (d)' DT50_c_in_sediment log(2)/kaas(5,3) log(2)/kaas(5,2)}];
tbl = [tbl;{'DT50 metabolite_2 s (d)' DT50_m2_in_sediment log(2)/kaas(6,3) log(2)/kaas(6,2)}];
tbl = [tbl;{'frac sub to m1' f_sub_to_m1_w kaas(7,2) kaas(7,3)}];
tbl = [tbl;{'frac sub to m2' f_sub_to_m2_s kaas(8,2) kaas(8,3)}];
tbl = [tbl;{'logKom sub' logKom [] []}];
tbl = [tbl;{'Diss w exp (d)' Half_max_water [] []}];
tbl = [tbl;{'sed/water ratio' om_to_water_ratio [] []}];
disp(TITLE);
disp(tbl);
This results in the following fit:
Schrap, S.M., G.L. Sleijpen, W. Seinen & A. Opperhuizen (1994). Sorption kinetics of chlorinated hydrophobic organic chemicals: Part I. The use of first-order kinetic multi-compartment models, Environmental Science And Pollution Research International (1), pp. 21-28
©2004 Henk Verhaar
back to main site