Proposed guidance on how aged sorption studies for

Proposed guidance on how aged sorption studies for
pesticides should be conducted, analysed and used in
regulatory assessments
S Beulke, W van Beinum
The Food and Environment Research Agency
Sand Hutton, York, YO41 1LZ, UK
J Boesten, M ter Horst
Alterra
PO Box 47, 6700 AA Wageningen, The Netherlands
April 2010
Funded by DEFRA within project PS2235
Table of contents
Table of contents ............................................................................................................................................... 2
1 Introduction ................................................................................................................................................... 3
2 Conceptual definition of equilibrium sorption ............................................................................................... 3
3 Experiments to derive aged sorption parameters ........................................................................................ 4
3.1
Soil selection and preparation ........................................................................................................... 4
3.2
Sample preparation and incubation ................................................................................................... 4
3.3
Sampling time points ......................................................................................................................... 5
3.4
Extraction and analysis ...................................................................................................................... 6
4 Fitting of kinetic models to data from aged sorption studies ........................................................................ 6
4.1
Data issues ........................................................................................................................................ 6
4.1.1 Data requirements ................................................................................................................. 6
4.1.2 Data handling ........................................................................................................................ 7
4.1.3 Outliers .................................................................................................................................. 7
4.2
Models ............................................................................................................................................... 8
4.3
Tools ................................................................................................................................................ 10
4.3.1 PEARLNEQ ......................................................................................................................... 11
4.3.2 ModelMaker 4.0 ................................................................................................................... 12
4.3.3 MatLab ................................................................................................................................. 12
4.4
Optimisation procedure.................................................................................................................... 12
4.4.1 Variables used in the optimisation....................................................................................... 12
4.4.2 Fitted parameters ................................................................................................................ 13
4.4.3 Optimisation settings ........................................................................................................... 13
4.4.4 Starting values ..................................................................................................................... 15
4.4.5 Parameter ranges ................................................................................................................ 16
4.4.6 Weighting ............................................................................................................................. 16
4.5
Goodness of fit and acceptance criteria .......................................................................................... 18
4.5.1 Visual assessment............................................................................................................... 18
2
4.5.2 Chi -test ............................................................................................................................... 20
4.5.3 Confidence intervals and relative standard error ................................................................ 21
4.5.4 Additional acceptance criteria ............................................................................................. 21
5 Use of aged sorption parameters in regulatory exposure assessments .................................................... 21
5.1
Sources of input data for regulatory exposure assessments .......................................................... 21
5.1.1 Estimation of sorption parameters from soil or pesticide properties ................................... 21
5.1.2 Default values ...................................................................................................................... 22
5.1.3 Experimental laboratory incubation studies ........................................................................ 22
5.1.4 Column and field studies ..................................................................................................... 22
5.2
Aged sorption in the tiered pesticide leaching assessment ............................................................. 22
5.3
Special considerations for metabolites ............................................................................................ 24
6 References ................................................................................................................................................. 25
Appendix 1: Fitting of a two-site aged sorption model with PEARLNEQ to two example datasets ................ 26
Draft - April 2010
Page 2 of 40
1
Introduction
Sorption of a pesticide to soil constituents determines its availability to non-target organisms and its potential to
move to groundwater or surface waters. Pesticide fate modelling at the first tier of regulatory environmental risk
assessments assumes that pesticide sorption is instantaneous and fully reversible. This implies that sorption
coefficients are constant with time. However, sorption has frequently been observed to increase with
increasing time of interaction with the soil (e.g. Walker and Jurado-Exposito, 1998; Cox and Walker, 1999).
Research for DEFRA project PS2206 (DEFRA, 2004) and PS2228 (DEFRA, 2009) confirmed that amounts
of pesticide in the soil solution are constantly changing.
Experimental studies that demonstrate an increase in pesticide sorption with time are increasingly submitted
to regulatory authorities as part of the regulatory data package. The results of these studies are used by
applicants to revise estimates of predicted environmental concentrations in groundwater. Pesticide leaching
models that include changes in sorption with time are used for this purpose. There is currently a lack of
agreed and clear guidance on how aged sorption studies should be conducted, analysed, interpreted and
hence used in regulatory exposure assessments. This document addresses this need. The guidance
outlined below was developed within a project sponsored by the UK Chemicals Regulation Directorate
(DEFRA, 2010). It will be presented and discussed at a workshop with representatives of European
regulatory authorities, academia, consultancies and industry. Revisions will be made based on the feedback
received.
The proposed guidance given in this document is underpinned by a literature review, experimental work and
modelling undertaken within the research project. Details are given in the research report (DEFRA, 2010). It
is important to note that the draft guidance is based on current knowledge. It is expected that some of the
recommendations will be modified in the future as new evidence comes to light. The need for revision of
certain aspects may become apparent when this draft guidance is put into practice.
2
Conceptual definition of equilibrium sorption
Sorption kinetics of pesticides in soils takes place at different time scales. Wauchope et al. (2002) distinguish
three time scales: (i) minutes, (ii) hours and (iii) weeks or years. Sorption increases very rapidly during the
first days after application. This is followed by a more gradual increase in sorption over time. Sorption over
the whole timescale can only be described accurately with multi-site models that conceptualise several types
of non-equilibrium sites reacting at different rates. These models have a large number of parameters and
more simplified two-site models with an equilibrium site and a single non-equilibrium site are preferred within
the regulatory context. Two-site models can either describe the initial rapid increase in sorption over the first
hours and days or the longer-term behaviour of sorption. But it is often difficult match sorption over the whole
timescale with a two-site model. Pesticide movement to depth by chromatographic leaching is mainly driven
by the sorption behaviour of the pesticide over the time scale of days to months. A two-site model that can
describe the increase in sorption from a few days after application onwards was therefore considered best
for regulatory leaching modelling.
A definition of the equilibrium fraction of the two-site model needs to be made for operational reasons. In the
model, the defined equilibrium fraction determines the initial sorption immediately after application. In this
guidance, the equilibrium fraction is defined as sorption measured during shaking of the soil with aqueous
solution for 24-hours. Sorption in soil at natural moisture conditions is initially lower than that estimated from
shaken 24-hour batch experiments. It may take approximately one week before the 24-hour value is
reached. However, sorption during the first week is expected to be less important for leaching to
groundwater than long-term sorption. Therefore it is probably justified to assume that the initial sorption
equals the amount of sorption in a 24-h shaken batch experiment. The operational definition recommended
here was also adopted by the current FOCUS groundwater scenarios work group (FOCUS, 2010). It is
consistent with the general perception that sorption equilibrium is reached within 24-48 hours.
The use of the 24 hour batch value as an operational definition of equilibrium sorption is more appropriate for
the description of pesticide losses to groundwater than to surface water. Entry into surface waters via
drainflow or runoff is mainly determined by short-term response to rainfall soon after application of pesticides
and less affected by long-term sorption. This is particularly true where preferential flow is an important
process. In this case, movement to drains can occur within the first hours or days of application and a correct
description of sorption at this time is important.
Draft - April 2010
Page 3 of 40
3
Experiments to derive aged sorption parameters
A standardised protocol to measure time-dependent sorption parameters for regulatory use must ensure the
reproducibility of the experimental results and maximise the reliability of derived model parameters. The
selection of the recommended procedure was based on a review of methods and experimental work
described by DEFRA (2010). A laboratory method was chosen because it is a well-defined system and
provides consistent and repeatable results that are relatively easy to interpret. Column and field studies are
more similar to the conditions of pesticide use in practice, but the greater complexity of these studies leads
to larger uncertainties in the model parameters and hence in the regulatory risk assessment. These methods
are thus not recommended. For a more detailed discussion see DEFRA (2010).
In brief, the recommended method is a laboratory incubation study where soil samples are treated with the
test substance and incubated in the dark at constant temperature and soil moisture. After selected time
intervals, samples are extracted with aqueous solution (to determine the concentration in the liquid phase
and extracted with solvent to determine the total extractable residue in the samples. The procedure
described below is similar to that recommended by OECD guideline 307 for aerobic and anaerobic
transformation in soil (OECD, 2002) except that an aqueous extraction step is added for measuring
desorption. A standard adsorption test (OECD 106, 2000) should be performed on the same soil to derive
the equilibrium sorption parameters.
To avoid duplication of effort, it is suggested that the additional measurements for aged sorption could be
routinely included in standard rate of degradation studies (OECD 307). The measurements would then be
available for modelling at the higher tier if required. To avoid the need for additional batch sorption studies, it
is recommended to use the soils selected for the standard OECD 106 batch sorption tests in the
degradation/aged sorption experiments. Instead of initiating aged sorption studies when the need for these
experiments becomes apparent in the lower tier risk assessment, it is proposed to include aged sorption
measurements in the routine suite of regulatory fate studies from the outset. Although this procedure will in
some cases generate work that will prove unnecessary, it will save considerable time and effort in those
cases where information on aged sorption is required.
3.1
Soil selection and preparation
• It is difficult to recommend a minimum number of aged sorption studies that must be undertaken. The
large variability in parameters from studies with the same pesticide applied to different soils and the
strong sensitivity of leaching models for aged sorption parameters suggests that the number of studies
should be large (>10). However, the experimental and modelling effort is substantial. It is thus
recommended to carry out a minimum of four aged sorption studies with contrasting soils.
Batch sorption is usually measured in five soils according to the guidance in OECD 106 (OECD, 2000).
The route and rate of degradation is measured in one soil and the rate of degradation is measured in
three additional soils as described in OECD guideline 307 (OECD, 2002). As there are no detailed
specifications of the soil properties for the three additional soils in OECD 307, it should be possible to use
the same soils in the degradation / aged sorption studies as in the batch sorption studies.
• Soils for the incubation study should be freshly sampled from the field and gently dried to a suitable
moisture content for sieving. The soils should be sieved through a 2-mm mesh. Soils may be stored for a
maximum of 3 months before the start of the incubation study, in a dark and cool place while protected
from drying out and leaving an opening to ventilate (i.e. same procedure as in OECD 307).
• Soil properties that need to be determined include the organic matter content, texture and pH of the soil
(see OECD 307). Preferably, the water retention at pF2 to 2.5 should be measured on a sieved soil
sample. Alternatively, the maximum water holding capacity of the sieved soil can be determined. The
actual moisture content of the soil is measured by oven-drying of subsamples just before the start of the
incubation study.
3.2
Sample preparation and incubation
• The soil is acclimatised at the appropriate moisture content (just below the target moisture content) and
temperature for at least 2 days, but ideally 7 days before applying the test compound. This can be either
before or after weighing out the individual soil samples. The soil is incubated in the dark at a constant
temperature of 20°C (± 2 °C). The moisture content should be between pF 2 and 2.5. Alternatively, the
water content can be adjusted to between 40 and 60% of the maximum water holding capacity (mwhc).
The moisture should be chosen such that the impact on the structure and aeration during mixing of the
Draft - April 2010
Page 4 of 40
soil with the pesticide is minimised. The moisture content should be sufficient to ensure adequate
microbial activity. It should be noted that there is no simple relationship between water retention and
maximum water holding capacity. The two ranges given here (pF 2 and 2.5 and 40-60% mwhc) are thus
not equivalent.
• The moisture content of the soils during sample preparation and incubation is maintained by monitoring
weight losses and by adding water when necessary.
• Individual samples are prepared for each sampling time. At least 2 replicates should be prepared for each
sampling time. More replicates may be needed if a large variation between replicates is expected.
Otherwise it may be more efficient to increase the number of sampling time points instead.
There are two general procedures to generate replicate samples. Both procedures are suitable to
generate true, independent replicates:
a. A large soil sample is treated with the pesticide, mixed and sub-samples are weighed into individual
flasks for incubation;
b. Smaller soil samples are treated individually with the pesticide and incubated in separate flasks.
• About 50 to 200 g of soil is weighed into each replicate incubation flask or jar.
• The application solution may be prepared using analytical grade chemical or formulated product,
whichever is appropriate. The application solution of the pesticide may be prepared in water or with small
amounts of solvent as described in the OECD guideline 307 (OECD, 2002). Only a minimum of solvents
should be used if necessary.
• Radiolabelled chemical may be used for the experiments as long as the percentage of parent compound
in the sample extracts is determined to calculate the residue and concentration of the parent compound.
• The selected application rate should correspond to the recommended application rate in the field,
assuming that in the field the applied compound would be mixed into the top 2.5 cm soil (unless
incorporated to a larger depth) as described in the OECD guideline 307 for degradation (OECD 2002).
• The application solution is distributed drop-wise on the soil surface to get an even distribution and to
maximise contact with the soil. The soil is left for one or two hours to absorb the liquid. Then the soil
samples are mixed gently with a spatula, taking care not to disrupt the soil structure as far as possible.
The weights of the flasks are noted so that the moisture content of the soils can be monitored.
• Alternatively the pesticide solution may be applied and mixed into a larger amount of soil (e.g. 1 or 2 kg)
and separated over several soil samples.
• The incubation period should normally not exceed 120 days as specified in the guidance for degradation
studies (OECD 307).
3.3
Sampling time points
At the selected time points, replicate samples are removed from the incubator and sacrificed for aqueous
and solvent extraction.
• Time intervals should be chosen so that the pattern of decline of the mass and aqueous concentration of
the test substance can be established. Time points should be closer together at the beginning of the
experiment and further apart towards the end of the experiment. At least six time points are needed for
the derivation of time-dependent sorption parameters. It should be noted that some time points may be
eliminated during the analysis of the raw data (Section 4). Six time points must remain thereafter.
• The earliest time point that should be included in the time-dependent sorption study is 48 hours. Earlier
time points are not recommended for the modelling (see Section 4). At least one sampling time point
should be within 3 days after application.
• For compliance with the guidelines for degradation studies (OECD 307) one additional sampling should
be undertaken immediately after application (0-day sample). This measurement should not be included in
the modelling for time-dependent sorption (see Section 4) therefore it is not necessary to extract the
sample with aqueous solution but it is sufficient to extract the sample with solvent only.
Draft - April 2010
Page 5 of 40
3.4
Extraction and analysis
The aqueous extraction is performed by gently shaking the soil with a solution of CaCl2 (0.01M) for 24 hours.
Then the samples are centrifuged and the concentration of parent compound is analysed in the supernatant.
The soil is extracted with solvent to determine the total extractable residues of the parent compound.
Aqueous extraction and solvent extraction may be performed consecutively on the same sample or in
parallel on sub-samples from the same flask. It is not appropriate to measure total and aqueous extractable
residues in samples that have been dosed separately.
The aqueous phase concentration must be characterised by shaking with CaCl2 for 24 hours. It is not
permitted to extract the soil water held by the moist soil during incubation by centrifugation. For a justification
of this recommendation, see DEFRA (2010).
• The soil samples need to be mixed well with a spatula before sub-samples are taken from the flasks. If
parallel samples are used for aqueous and solvent extraction then both sub-samples need to be taken
from the same flask.
• For the aqueous extraction, the soil is extracted by shaking with CaCl2 solution (0.01M). The soil:solution
ratio should be chosen based on the soil:solution ratio in the batch sorption experiment on the same soil .
The soil is shaken gently for 24 hours at the lowest rate possible at which the soil would stay suspended
in the liquid and no solids are settling on the bottom of the tube. The low speed is required to keep the
disruption of the soil structure during aqueous extraction to a minimum. Then most of the supernatant is
removed after centrifugation and the concentration of parent compound is analysed in the liquid.
• Then the samples are extracted with solvent to determine the extractable residues of the parent
compound. The solvent extraction method must recover 90-110% of the applied compound just after
application. This range applies to radiolabelled and non-radiolabelled studies. A larger deviation cannot
be permitted because this leads to errors in the estimated model parameters. The same method should
be used throughout the experiment irrespective of the extraction efficiencies at later time points. It is
important to reflect on which fraction of the pesticide should be captured by the solvent extraction in an
aged sorption study. Harsher extraction methods will lead to larger aged sorption parameters (i.e. less
conservative with respect to leaching) but slower (i.e. more conservative) degradation parameters.
• The concentration of the parent compound in the aqueous extract and the total extracted mass of parent
compound in the soil are should be determined. If consecutive extraction is used then both extracts need
to be accounted for in the calculation of the total extractable residue.
• The limit of quantification (LOQ) for the parent compound should be determined in aqueous and solvent
extracts. Measurements below the LOQ are not included in the modelling (see Section 4).
4
Fitting of kinetic models to data from aged sorption studies
4.1
Data issues
4.1.1
Data requirements
The quality of the dataset and the handling of the data influence the estimated sorption parameters. The
following minimum requirements should be met:
• The incubation study should follow the guidance given by OECD 307 (OECD, 2002) and the additional
recommendations given in Section 3 of this guidance document. Batch sorption studies to determine the
Freundlich exponent N must be undertaken in accordance with OECD 106 (OECD, 2000)
• The system must be well characterised. The mass and water content of the soil during incubation, the
volume of water added during extraction, the duration and intensity of the extraction should be stated.
Information on the texture, organic carbon content, pH and water retention or maximum water holding
capacity of the sieved soil should also be available.
• Data on total mass and aqueous concentration must be available. The total mass is defined as the mass
that is extractable by organic solvents. The model considers non-extractable residues to be equivalent to
transformation products, the non-equilibrium sorption component is independent of the mechanism by
which the compound is ‘lost’ from the system. Measurements of solvent-extractable pesticide in % of
applied radioactivity are suitable if the radioactivity is characterised.
• Experimental studies must provide sufficient and adequate sampling points to ensure a robust estimation
of parameters. The pattern of decline in mass and concentration must be well established.
Draft - April 2010
Page 6 of 40
• It is often difficult to obtain a good fit of the model to the data over the whole timescale of the experiment.
Measurements within the first two days after application are often not matched by the model. Processes
other than long-term sorption, such as short-term adsorption, precipitation and dissolution are likely to
influence the measurements of both mass and concentration during the first two days after pesticide
application. The two-site model is not able to describe kinetically the rapid reactions that occur within the
first hours and days after application as well as the slower processes operating at a time scale of weeks
or months because it contains only one kinetic sorption site (see also Section 2). In order to describe the
long-term behaviour well, samples that were taken less than 48 hours after treatment should be excluded
from model fitting.
• The first sampling time that is included in the modelling must be between 48 and 72 hours after treatment.
Sorption increases rapidly within the first days of the experiment. Initial sorption is estimated by the
model. An accurate estimate cannot be obtained if the first measurement is taken later than 72 hours
after treatment. The fitting may not be reliable for rapidly degrading compounds where the total mass
declines quickly between application and the first acceptable measurement. The fitted initial mass should
thus always be compared with the amount measured immediately after application. A discrepancy of
more than 15% is not acceptable as this will introduce a bias into the other model parameters. If the mass
at the time of treatment is not measured, then the fitted initial mass can be compared with the added
amount. It should, however, be noted that the discrepancy is then also dependent on the experimental
recovery inherent in the extraction method.
4.1.2
Data handling
• The measurements in the aged sorption study and the batch sorption study must not be corrected for the
recovery of the test compound.
• Measured data should be reported with a precision of at least 4 significant figures.
• The time of sampling must be rounded to one digit after the decimal point (e.g. 5.1 days, not 5.125 days).
• Incubation studies should be carried out with at least two true, independent replicates. Replicate values
for each sampling interval should not be averaged before curve fitting. Replicate analytical results from a
single sample are not true, independent replicates and should be averaged and treated as one sample
during parameter optimisation.
• Experimental results often include measurements below the limit of quantification (LOQ). Measurements
below the limit of quantification (LOQ) are uncertain and these should be discarded. If one of the replicate
measurements is missing or discarded because the value is below LOQ, then all measurements on this
sampling date and all subsequent dates must be discarded for both mass and concentration. This
deviation from guidance by FOCUS (2006) is necessary because the measurements are weighted during
the model fitting (see Section 4.4.6). The weight is equal to 1/measurement. This gives small
measurements a very large weight and these have a critical influence on the fitted aged sorption
parameters. Values below LOQ are not determined with sufficient precision and these must be excluded
from the fitting.
• A robust optimisation of parameters is only possible if the number of observations is appreciably larger
than the number of model parameters. The total number of sampling dates remaining after the elimination
of measurements <48 hours, measurements below the limit of quantification or outliers must not be
smaller than six.
• It is recommended to calculate the sorbed amount of pesticide for each sampling time as the difference
between the total mass and the mass in the liquid phase. Apparent linear sorption coefficients (Kd app)
can then be calculated as the ratio of sorbed:dissolved concentration. Kd app values will not be used in
the optimisation, but this variable can be useful in the interpretation of the data (see Section 4.5.1).
4.1.3
Outliers
Outliers in laboratory studies can be individual or several replicates or sampling dates. Outliers that are
explained by experimental errors should be eliminated before curve fitting.
Measurements that strongly differ from others without any obvious experimental reason should initially be
included in the optimisation. They can then be eliminated based on expert judgement and the fitting
procedure can be repeated. The results for the fits with and without outliers must be reported. The removal
of any data points as outliers must be clearly documented and justified in the report.
Draft - April 2010
Page 7 of 40
If one replicate measurement for mass and/or concentration is an outlier, only this needs to be eliminated
and the other replicate(s) can be included in the optimisation. After elimination of outliers, there may be more
measurements for the total mass than the aqueous concentration or vice versa. The discrepancy must be
small to ensure that the fit is not influenced much more strongly by one variable than the other.
4.2
Models
A number of models exist to describe aged sorption of pesticides in soils (see DEFRA, 2010, Section 2).
Only two-site adsorption-desorption models are currently considered suitable for regulatory use because
they provide a reasonable balance between the complexity of the model and the experimental effort required
to determine the model parameters. Two-site models are now implemented into the software packages
FOCUS PEARL, MACRO 5.0 onwards, FOCUS PELMO and FOCUS PRZM to enable the simulation of
1
kinetic sorption (FOCUS, 2010).
FOCUS PEARL
The leaching model FOCUS PEARL uses the two-site model according to Leistra et al. (2001). The same
two-site model is implemented for a laboratory system in the PEARLNEQ software. This software can be
used to derive input parameters for FOCUS PEARL. The PEARLNEQ model is depicted in Figure 4-1.
Figure 4-1. Schematic representation of the PEARLNEQ model showing the soil solution on the right and the
equilibrium and non-equilibrium sorption sites on the left. Only pesticide in the equilibrium domain (indicated by
the dashed line) is subject to degradation.
Freundlich:
KF,EQL nF
equilibrium
Freundlich
KF,NEQ nF
non-equilibrium
sorption
sorption
Desorption Rate Constant:
kdes
Ratio KF NEQ:KF EQL
The model assumes that sorption is instantaneous on one fraction of the sorption sites and slow on the
remaining fraction (Leistra et al., 2001). The model does not account for irreversible sorption. Degradation is
described by first-order kinetics. Only molecules present in liquid phase and sorbed to the equilibrium site
are assumed to degrade. Molecules sorbed on the slow non-equilibrium sorption site are considered not to
degrade. The PEARLNEQ model can be described as follows:
M p = V c L + M s ( X EQ + X NE )
X EQ
 c 
= K F , EQ cL , R  L 
 cL , R 
(1)
N
(2)
N
 c 
dX NE
= kd ( K F , NE cL , R  L  − X NE )
dt
 cL , R 
(3)
K F , NE = f NE K F , EQ
(4)
dM p
dt
= −k t (V c L + M s X EQ )
KF,EQ = mOM KOM,EQ
(5)
(6)
1
Note that at the time of writing this document, the versions of FOCUS PELMO and FOCUS PRZM that include aged
sorption are not officially released.
Draft - April 2010
Page 8 of 40
where:
Mp
V
Ms
cL
cL,R
XEQ
XNE
KF,EQ
KF,NE
N
kd
fNE
kt
mOM
KOM,EQ
=
=
=
=
=
=
=
=
=
=
=
=
total mass of pesticide in each jar (µg), acronym Mas
the volume of water in the soil incubated in each jar (mL), acronym VolLiq
the mass of dry soil incubated in each jar (g), acronym MasSol
concentration in the liquid phase (µg/mL), acronym ConLiq
reference concentration in the liquid phase (µg/mL), acronym ConLiqRef
content sorbed at equilibrium sites (µg/g)
content sorbed at non-equilibrium sites (µg/g)
equilibrium Freundlich sorption coefficient (mL/g), acronym CofFreEql
non-equilibrium Freundlich sorption coefficient (mL/g), acronym CofFreNeq
Freundlich exponent (-), acronym ExpFre
-1
desorption rate coefficient (d ), acronym CofRatDes
a factor for describing the ratio between the equilibrium and non-equilibrium Freundlich
coefficients (-), acronym FacSorNeqEql
-1
= degradation rate coefficient (d )
= mass fraction of organic matter in the soil (kg/kg), acronym CntOm
= coefficient of equilibrium sorption on organic matter (mL/g), acronym KomEql
The model has six parameters: the initial concentration of the pesticide, the degradation rate constant kt, the
equilibrium sorption coefficient KOM,EQ, the Freundlich exponent N, the ratio of non-equilibrium sorption to
equilibrium sorption fNE and the desorption rate constant kd. Note that the notation N for the Freundlich
1
exponent used here is equivalent to the /n that is commonly used in the Freundlich equation.
MACRO
A very similar model has been implemented into the pesticide leaching model MACRO (Larsbo and Jarvis,
2003). It is based on the model by Streck et al. (1995). The rate equation used by PEARLNEQ (Equation 3)
differs from that used by MACRO:
N
 c 
dX NE
α
= MACRO ( K F ,Total cL, R  L  − X NE )
dt
f NE MACRO
 cL, R 
(7)
The definition of fNE is also different in MACRO. Here, fNE expresses non-equilibrium sorption as a fraction of
total sorption (Equation 8) whereas fNE in PEARLNEQ is the ratio of non-equilibrium to equilibrium sorption
(Equation 4).
f NE MACRO =
K F , NE
K F , EQ + K F , NE
(8)
where:
XNE
αMACRO
fNE MACRO
KF,Total
KF,EQ
KF,NE
=
=
=
=
=
=
content sorbed at non-equilibrium sites (µg/g)
-1
desorption rate coefficient (d ) used in MACRO.
fraction of the non-equilibrium sorption sites in MACRO (-)
sum of equilibrium plus non-equilibrium Freundlich sorption coefficient (mL/g)
equilibrium Freundlich sorption coefficient (mL/g)
non-equilibrium Freundlich sorption coefficient (mL/g)
The degradation rate on the non-equilibrium sites in MACRO can be set equal to the rate in the equilibrium
pool or to zero. Zero degradation in the non-equilibrium pool is identical to the concepts in PEARLNEQ. The
relationship between the parameters used in MACRO and PEARLNEQ (FOCUS, 2010) is:
f NE MACRO =
f NE PEARL
1 + f NE PEARL
α MACRO = kd PEARL
Draft - April 2010
and
f NE PEARL
1 + f NE PEARL
f NE PEARL =
and
f NE MACRO
1 − f NE MACRO
k d PEARL =
α MACRO
f NE MACRO
(9, 10)
(11,12)
Page 9 of 40
PELMO and PRZM
The FOCUS GW II versions of PELMO and PRZM use the same model as FOCUS PEARL. The parameters
of the PEARLNEQ or Streck model can be entered into PELMO. The Streck values are then converted within
the model.
4.3
Tools
Several tools are available for fitting the two-site model to the data. The model parameters are derived by an
optimisation procedure. The estimation of parameter values from aged sorption studies consists of several
steps:
1. Entering the measured data for each sampling time.
2. Making an initial guess for each parameter value of the selected model (referred to as “starting value”).
3.
Calculation of the data at each time point.
4.
Comparison between the calculated and measured data.
5. Adjustment of the parameter values until the discrepancy between the calculated and measured
concentrations is minimised (“best fit”).
Steps 3-5 are carried out automatically within software tools. These packages start from the initial guess
made by the modeller and repeatedly change the parameter values in order to find the best-fit combination.
In order to use such an automated procedure, “best fit” must be defined in the form of a mathematical
expression referred to as ‘objective function’. Often, the sum of the squared differences between the
calculated and observed data (sum of squared residuals = SSQ) is used. The software package aims at
finding the combination of parameters that gives the smallest SSQ. This method is referred to as least
squares method. Maximum likelihood methods can also be used. These maximise the probability that the
simulated curve is an exact match of the measured data.
The method to adjust the parameter values from the previous guess based on the objective function differs
between different tools. Many optimisation packages use the Levenberg-Marquardt algorithm. This method
linearises the differential model equations and calculates the model output for the initial parameter guess
based on the linear equation. It then changes the parameters one at a time up or down (or in both
directions), calculates the model output again and compares the objective function between the old and new
parameter value(s). The change in the objective function drives the size and direction of the next change in
the parameter value. When the objective function does no longer change, the parameter value at that point is
returned as the optimum value. The standard error of the parameter is calculated as a function of i) the value
of the objective function at the optimum, ii) the total number of observations, iii) the number of parameters
and iv) the linearised form of the differential equations. The confidence interval is calculated from the
standard error based on the assumption that the standard errors are normally distributed.
An alternative approach is the Markov Chain Monte Carlo method. The Levenberg-Marquardt algorithm
varies parameters within the constraints specified by the user and gives equal probability to all values
between these boundaries. In contrast, the expected type and width of the parameter distribution can be
specified in the Markov Chain Monte Carlo method. For example, it may be expected the parameter DegT50
lies somewhere within a log-normal distribution with a mean of 20 days and a standard deviation of 5. This
gives values near 20 a higher probability than values at the tails of the distribution. A parameter value is
selected from this distribution and the objective function is calculated. The parameter value is then changed
and the objective function is calculated again. The parameter distribution is updated during the optimisation
based on the differences between the objective functions at each step. The final distribution gives
information on the most likely parameter value that gives the best fit. The confidence intervals can be derived
directly from the final parameter distribution.
The Levenberg-Marquardt algorithm changes the parameter value up or down from its starting point. It can
get ‘trapped’ in a region where the objective function is small (’local minimum’) without realising that even
smaller objective functions (‘global minimum’) could be achieved if the parameter jumped to a value far away
from the starting point. The Markov Chain Monte Carlo method evaluates the objective function for the whole
distribution of possible parameter values. It is, thus, in principle more likely to find the global minimum of the
optimisation than the Levenberg-Marquardt algorithm, provided the prior distribution includes the true
optimum parameter. However, the settings for the Levenberg-Marquardt algorithm can be fine-tuned to
ensure that the global minimum is reached.
Draft - April 2010
Page 10 of 40
Three tools that are commonly used to derive aged sorption parameters are briefly described below.
Alternative optimisation packages can be used provided the tool and optimisation settings give robust fits.
The independence of the optimised parameter values from the starting values must be demonstrated
because this increases the likelihood that the global minimum can be reached. The optimisation package
must also provide the output that is required to assess the goodness of fit according to Section 4.5 (e.g.
confidence interval or standard error). Ideally, the results from the alternative tool should be compared with
those from one of the three tools described below. This is intended to be a one-off test of the alternative
optimisation package, a comparison with other tools is not required after the similarity of results has been
demonstrated for example datasets.
4.3.1
PEARLNEQ
PEARLNEQ combines the two-site sorption model that is implemented in FOCUS PEARL with the
optimisation software PEST (Doherty, 2005). The model is simultaneously fitted against data on the total
mass of the pesticide in soil (µg) and the concentration in the liquid phase (µg/mL). PEARLNEQ is run
repeatedly by PEST and the parameters are adjusted until the best possible fit to the measured data is
achieved based on the least squares method and the Gauss-Marquardt-Levenberg algorithm. The program
is DOS based and operates on command file or command line level. Boesten et al. (2007) provide a short
description of PEARLNEQ.
The program package of PEARLNEQ includes the PEARLMK.EXE program that produces all necessary
PEST files with the help of a text file with the extension .mkn. In order to carry out the non-equilibrium
parameter estimation procedure in PEARLNEQ, the *.mkn file of the PEARLNEQ package has to be
compiled following the instructions in the PEARLNEQ manual. The *.mkn file of PEARLNEQ for an example
case is given in Appendix 1.
The output generated by PEST includes the fitted parameters and their 95% confidence intervals, the sum of
squared residuals and daily output of the calculated total mass and liquid phase concentration for a period
specified by the user.
Some changes to the PEST input files and the PEARLNEQ code have been implemented into the latest
version of PEARLNEQ (Version 5) to comply with the guidance given in this document. The following
2
modifications were made in Version 5 :
2
•
PEARLNEQ Version 4 optimises the parameters Mp ini, kt, kd and fNE. The equilibrium Freundlich
sorption coefficient KOM,EQ must be specified by the user and is not included in the optimisation. This
parameter is usually set to sorption coefficient from a standard 24-hour shaken batch study. DEFRA
(2010) found that the 24-hour batch sorption coefficient is sometimes a poor estimate of sorption at
the beginning of the incubation study. PEARLNEQ Version 5 was modified to enable users to
optimise KOM,EQ against the measured data.
•
The optimisation settings in PEST were found to result in variable results depending on the starting
values. The settings were changed to those specified in Section 4.4.3.
•
The time step for the numerical solution of the differential equations is set to a constant value by the
user in Version 4. The time-step varies automatically depending on the rate constants kt and kd in
Version 5.
•
The PEST files were adjusted such that replicate measurements can be fitted simultaneously.
•
In Version 4, the measurements can only be entered for whole days. Fractions of a day to a
precision of one digit after the decimal point (e.g. 5.1 days) will be accepted by Version 5.
•
Version 4 provides daily model output for the total mass, liquid phase concentration and the sorbed
mass in the non-equilibrium phase. Output for the sorbed mass in the equilibrium phase was added
in Version 5. Output for the apparent linear sorption coefficient Kd app (ratio of sorbed : dissolved
concentration) is also given in Version 5. This variable will not be used in the fitting procedure, but it
can be useful in the assessment of the visual fit. All output is now generated for each tenth of a day.
At the time of writing this document, Version 5 is not publicly available
Draft - April 2010
Page 11 of 40
4.3.2
ModelMaker 4.0
TM
ModelMaker is one of the tools that are recommended for parameter fitting within the framework of FOCUS
kinetics (a more detailed description can be found in FOCUS, 2006). It allows users to build their own
models using inter-linked variables or compartments. Gurney and Hayes (2007) describe an implementation
TM
TM
(Figure 4-2). ModelMaker
of the two-site aged sorption model by Leistra et al. (2001) into ModelMaker
allows the user to optimise the equilibrium sorption coefficient KOM,EQ. Several replicates can be fitted
simultaneously.
TM
ModelMaker provides output of the optimised parameter values and their standard error, a graphical plot of
the measured and calculated data and the calculated values in tabulated form.
The ModelMaker
TM
version of the two-site aged sorption model can be made available on request.
TM
Figure 4-2. Implementation of non-equilibrium sorption in ModelMaker
4.3.3
MatLab
TM
MatLab (2007) is a numerical computing environment and fourth generation programming language.
TM
Developed by The MathWorks®, MatLab
allows matrix manipulation, plotting of functions and data,
implementation of algorithms, creation of user interfaces, and interfacing with programs in other languages.
TM
MatLab can be applied to build and solve mathematical models such as the two-site sorption model. Addon toolboxes are available for solving differential equations and to solve the optimisation of model
TM
parameters. The MatLab code can be tailored to the user’s requirements.
TM
BayerCrop Science integrated the two-site sorption model into an Excel® spreadsheet that calls MatLab
via Excel Link™. The parameters are adjusted based on the least squares method and the MarquardtTM
Levenberg algorithm. This is an option within the MatLab routine lsqnonlin (Solve nonlinear least-squares
data-fitting problems ). The default optimisation settings are used. The Markov Chain Monte Carlo method
could be implemented instead of the Marquardt-Levenberg algorithm. Further modifications could be made
to bring the version in line with the guidance outlined in this document (e.g. fitting of KOM,EQ, additional
graphical outputs). The tool generates various statistical outputs.
The current FOCUS GWII group fitted the two-site aged sorption model to the total mass and liquid phase
TM
concentration for an example dataset using the three software tools PEARLNEQ, ModelMaker
and
TM
MatLab . The results for all three tools were almost identical (FOCUS, 2010).
4.4
Optimisation procedure
This guidance below refers to the optimisation of the aged sorption model by Leistra et al. (2001). The
procedures for the optimisation of the Streck model are very similar.
4.4.1
Variables used in the optimisation.
The two-site aged sorption model comprises several variables (total mass, mass sorbed in equilibrium
phase, mass sorbed in non-equilibrium phase, concentration in liquid phase). The model should ideally be
fitted to the data on total mass and concentrations in the liquid phase because these are directly measured
during the experiment. An alternative procedure was tested by the FOCUS GW II group (FOCUS, 2010).
MatLab was used to fit the two-site model to the sorbed mass in the equilibrium and equilibrium phase.
These variables were calculated from the measured organic solvent and aqueous extractable residues. The
parameters derived with this method were compared with those optimised against the total mass and
concentrations in the liquid phase. The FOCUS GW II group found that the parameter values were
Draft - April 2010
Page 12 of 40
independent of the variables fitted, but the standard deviation of the parameters was smaller for the fits to
sorbed mass. However, additional modelling showed that the two methods are equivalent.
In radiolabelled studies, the radioactivity measured in the aqueous and solvent extracts must be
characterised and converted to mass and concentrations of the parent compound or metabolite of interest.
4.4.2
Fitted parameters
Non-equilibrium sorption model
The two-site aged sorption model described by Leistra et al. (2001) has six parameters (Mp ini, KOM,EQ, N, kt,
kd and fNE), see Section 4.1. All parameters except N should be optimised against measured data. In the
optimisation tool PEARLNEQ, the parameter kt is not optimised directly. The degradation half-life (DegT50,
days) is optimised instead and kt is calculated within the model as ln(2)/ DegT50.
KOM,EQ should not be taken from the batch study because there is often a discrepancy between the
measurements in the batch experiment and the aged sorption study. This is due to (i) the variability inherent
in all experimental studies which leads to differences between repeat studies, and (ii) the small difference in
the metholodology (the pesticide is added to the CaCl2 shaking solution in the batch study whereas it is
dripped onto the soil or mixed into the soil prior to extraction in the aged sorption study).
Fixing KOM,EQ to the batch value could introduce a bias into the fitting of the aged sorption parameters. The
equilibrium Kom value in the aged sorption study can be estimated more accurately by curve fitting. The
fitted value should be compared with the batch value as a plausibility check. A discrepancy of more than
20% is not acceptable.
Ideally, the Freundlich exponent N should also be derived from the aged sorption study as it is not
guaranteed that N is the same in the batch and aged sorption experiment. However, a robust estimation of n
requires studies carried out at several initial pesticide concentrations. These are rarely available. The
optimisation of all six parameters against the data from an aged sorption study at a single concentration
would lead to large uncertainties in the results. The fit is considered to be more robust if N is taken from the
batch study and fixed during the optimisation. The batch study must be undertaken with the same soil as that
used in the aged sorption study. Values from studies with other soils (or averages of several soils) are not
acceptable.
Equilibrium sorption model
A model fit should also be undertaken with equilibrium sorption only. The non-equilibrium component of the
model can be switched off by fixing fNE and kd at zero. Only Mp ini, DegT50, and KOM,EQ are then optimised.
The results of this optimisation are used as a benchmark for comparison with the fit by the aged sorption
model.
4.4.3
Optimisation settings
The optimisation criterion (‘objective function’) is often the minimisation of the sum of squared residuals
between the measured data and the simulated values (SSQ). There may be a single combination of
parameters that results in the smallest possible value for the sum of squared residuals (“global minimum”).
But there are often several additional combinations that also result in small SSQs (“local minima”). In
particular, the parameters fNE and kd are strongly related. The increase in one of the two parameters can be
compensated to some extent by a decrease in the other parameter. Various combinations of fNE and kd may
thus result in similar fits. This is referred to as non-uniqueness. In this case, the software may stop the
optimisation procedure before the global minimum is found.
The ability to reach the global minimum depends on the initial guess (the closer the initial guess to the best
possible value, the better), the nature of the specific optimisation problem and the settings within the
software package. Different parameters may be obtained by different software packages and the derived
combination of parameters does not necessarily provide the best possible fit to the measured data.
The problem of non-uniqueness can be minimised by selecting certain optimisation settings. The
recommended settings in the PEST control file that is provided with the PEARLNEQ programme are given in
Table 4-1. For definitions of the PEST parameters see the user manual (Doherty, 2005).
Draft - April 2010
Page 13 of 40
Table 4-1. PEST control settings
PEST parameter
description
Value
PRECIS
Precision used when writing parameter values to model input files
(single or double)
single
DPOINT
Use of decimal point when writing parameter values to model input
files (point or nopoint)
point
RLAMBDA1
Initial lambda
5
RLAMFAC
Lambda adjustment factor
2
PHIRATSUF
Sufficient new/old phi ratio per optimisation iteration
0.1
PHIREDLAM
Limiting relative phi reduction between lambdas
1.0E-02
NUMLAM
Maximum trial lambdas per iteration
15
RELPARMAX
Maximum relative parameter change (relative-limited changes) (used
if PARCHLIM is ‘relative’)
na
FACPARMAX
Maximum factor parameter change (factor-limited changes) (used if
PARCHLIM is ‘factor’)
4
FACORIG
Fraction of initial parameter values used in computing; change limit for
near-zero parameters
1.0E-03
PHIREDSWH
Relative phi reduction below which to begin use of central derivatives
(used if FORCEN = ‘switch’)
na
NOPTMAX
Maximum number of optimisation iterations
50
PHIREDSTP
Relative phi reduction indicating convergence
0.10E-02
NPHISTP
Number of phi values required within this range
5
NPHINORED
Maximum number of consecutive failures to lower phi
10
RELPARSTP
Minimal relative parameter change indicating convergence
0.10E-02
NRELPAR
Number of consecutive iterations with minimal parameter change
4
INCTYP
Increment type (used if FORCEN = ‘always_2’ or ‘switch’)
na
DERINC
Increment (used if FORCEN = ‘always_2’ or ‘switch’)
na
DERNCLB
Increment lower bound (used if FORCEN = ‘always_2’ or ‘switch’)
na
FORCEN
Forward difference, central difference or both used in course of an
optimisation run (resp. always_2, always_3, switch)
always_3
DERINCMUL
Multiplier
2
DERMTHD
Variants of the central (i.e. three point) method of derivatives
calculation (‘parabolic’, ‘best_fit’, ‘outside_pts’)
best_fit
PARTRANS
Transformation (‘none’, ‘log’, ‘fixed’, ‘tied’)
none
PARCHGLIM
Change limit (‘relative’, ‘factor’)
factor
Draft - April 2010
Page 14 of 40
TM
The recommended optimisation settings in ModelMaker are shown in Figure 4-3. The accuracy of the
model integration (relative error per integration step) can be specified under Run Options (Model, Integrate,
-7
Advanced). It should be set to a small, very accurate value (e.g. 1×10 ).
Figure 4-3. Recommended optimisation settings in ModelMaker
TM
For other software tools please refer to the respective user manual.
4.4.4
Starting values
Different optimised values can be returned by the software for different combinations of initial guesses for the
parameters provided by the modeller (starting values). The optimisation settings specified above for PEST
TM
and ModelMaker
will reduce the dependency on starting values, but the problem of non-uniqueness can
not be fully overcome. The optimisation should thus be repeated with a number of different initial
combinations of parameter values. The results of all fits should be reported and the parameter combination
that gives the best objective function (e.g. the smallest SSQ) should be selected. If several starting values
give identical objective functions, then the combination with the smallest relative confidence intervals
(confidence interval as a fraction of the mean estimate) for fNE and kd should be chosen.
The following specific recommendations can be made:
• The initial mass Mp ini, is often close to the measured concentration at the first sampling point and
this can be used as a starting value in the optimisations where appropriate. An alternative is to use
the added mass. The starting value for the initial mass can also be derived by fitting a first-order
dissipation model to the data in a separate model run with any appropriate tool.
• The initial value for the degradation half-life DegT50 should be set to the first-order DT50 value. This
can be derived by fitting a first-order model to the data in a separate model run with any appropriate
tool.
• The batch KOM,EQ value measured in the same soil can be used as a starting value for the Freundlich
equilibrium sorption coefficient.
• At least four different initial guesses should be tested for fNE and kd (Table 4-2) The same starting
value for Mp ini, KOM,EQ and DegT50 can be used in all optimisations.
Table 4-2. Starting values for fNE and kd
fNE
0.2
0.2
1.5
1.5
Draft - April 2010
kd
0.004
0.05
0.004
0.05
Page 15 of 40
4.4.5
Parameter ranges
For some parameters, it may be useful to define ranges within which the parameter will be varied during
optimisation. This will prevent convergence at unrealistic local minima. A lower boundary > 0 will avoid
numerical problems during the optimisation (division by zero). It is recommended to constrain fNE between
0.001 and 10. The constraint range for kd should be 0.00001 and 0.5. Boundaries for Mp ini, KOM,EQ and
DegT50 may also need to be set.
4.4.6
Weighting
Aged sorption models should be simultaneously fitted to measurements for the total mass of a pesticide in
soil and the concentration in the liquid phase. The absolute values for the mass are often much larger than
the concentrations depending on the strength of sorption and the unit used (e.g. µg for the total residue and
µg/mL for the aqueous concentration). The same relative deviation of the modelled data from the calculated
values results in much greater squared residuals when the absolute value of the measurement is large. As a
result, an unweighted model fit will usually be dominated by the total mass and only marginally influenced by
the liquid phase concentrations. This can lead to a good fit to the mass, but a poor fit to the concentrations.
This can also result in large confidence intervals for the parameters kd and fNE.
The measurements must be weighted during the optimisation to minimise this problem. Weighted fitting
applies a correction factor to the residuals:
m
( )
Φ = ∑ wi ri
i =1
2
(13)
where Φ is the object function, ri is the residual (difference between the simulated and the measured value
corresponding to measurement i ), wi is the weighting factor and m is the total number of measurements
(sum of number of measurements of Mp and cL).
The preferred option is to define wi as the inverse of the measured value of Mp or cL. This will reduce the
weight of the mass data and increase the weight of the concentration data compared with unweighted fitting.
2
o ∆c

 ∆ M p, i 
 + ∑  L, j 
Φ = ∑




i =1  M p , i 
j =1  cL , j 
n
2
(14)
where n is the number of measurements for the mass and o is the number of measurements for the
concentration in the liquid phase (note that n=o unless outliers were eliminated), ∆Mp,i is the difference
between the simulated and observed mass for measurement i , Mp,i is the observed mass for measurement i,
∆cL,j is the difference between the simulated and observed concentration in the liquid phase for
measurement j, and cL,i is the observed concentration in the liquid phase for measurement j.
The time series of mass data consists of larger values at the beginning of the experiment and smaller values
at the end. The same is true for the time series of concentration data. Weighting by the reciprocal value
implies that the relative error in the measurements is constant with time, i.e. larger values for mass and
concentration are measured with the same relative accuracy than small values. This assumption was
supported by an analysis of measured data by DEFRA (2010). However, the number of datasets included in
the analysis was relatively small. If there is evidence that measurements at earlier time points are more
accurate than those at later time points (i.e. the relative error is smaller at earlier time points), then an
alternative weighting method could be used:
2
o ∆c
 ∆ M p, i 

 + ∑  L, j 
Φ = ∑


M p 
cL 
i =1 
j =1 
n
2
(15)
The weighting factor for the mass is here the inverse of the mass averaged over all time points M p whereas
the weighting factor for the concentration is the inverse of the concentrations averaged over all time points cL .
Both equation 14 and 15 account for the differences between mass and concentration measurements. The
difference is that equation 14 assumes the same relative error over time whereas equation 15 assumes that
the absolute error is constant over time. As a result, equation 14 gives a closer visual fit to the later time
points whereas equation 15 gives a closer visual fit to the earlier time points.
Draft - April 2010
Page 16 of 40
Weighting by 1/measurement (equation 14) is one of the options implemented in PEARLNEQ. The
TM
optimisation settings in ModelMaker
should be set to those shown in Figure 4-4 to match those in
PEARLNEQ (click on Advanced to access the weighting options).
TM
Figure 4-4. Recommended settings for data weighting in ModelMaker
– Option 1
Alternatively, the weights can be entered in the model data table as an additional column (Figure 4-5).
TM
ModelMaker divides the residuals by the weight specified by the user. The weights must thus be identical
to the measurements (and not the inverse value).
Figure 4-5. Recommended settings for data weighting in ModelMaker
Draft - April 2010
TM
– Option 2
Page 17 of 40
4.5
Goodness of fit and acceptance criteria
The decision on whether a model fit is acceptable or not should be based on:
• An assessment of the visual fit of a model with and without time-dependent sorption;
2
• A χ -test to assess the goodness of fit and to compare a model with and without time-dependent
sorption;
• An assessment of the confidence in the parameter estimates;
• The additional acceptance criteria specified in Section 4.5.4.
An example assessment of the goodness of fit is presented in Appendix 1.
4.5.1
Visual assessment
Measured and fitted data must always be presented graphically and a visual assessment of the goodness of
fit must be made. Model fits should be undertaken with and without non-equilibrium sorption. For each of the
two fits, four plots should be made (only the results for the starting values of fNE and kd that give the best fit
need to be plotted):
1.+2. Measured mass and aqueous concentration data and the calculated curves should be plotted versus
time.
3.+4. Calculated minus measured mass and calculated minus measured aqueous concentration data should
be plotted versus time (residuals). This is useful for revealing patterns of over- or under-predictions.
For an exact fit, all residuals are zero. Systematic deviations become apparent when negative and
positive residuals are not randomly scattered around the zero line.
Two additional plots can facilitate the interpretation of the results:
Apparent linear Kd values (Kd app) should be calculated from the measured data and the simulated
concentrations and plotted against time. Apparent Kd values are usually more scattered than the data on
mass and concentration. It is important that the apparent Kd value shows an increase over time that can be
distinguished from the scatter in the data. Figure 4-6 gives an example of an acceptable and unacceptable
pattern of Kd app.
35
16
30
14
25
12
Kd app (L/kg)
Kd app (L/kg)
Figure 4-6. Example of acceptable (left) and unacceptable (right) patterns of apparent Kd values
20
15
10
5
10
8
6
4
2
0
0
0
20
40
Time (d)
60
80
0
20
40
60
Time (d)
80
100
120
It is also important to compare the modelled line with the experimental apparent Kd values. Sometimes,
mass and liquid phase concentration are described well by the model, but the apparent Kd is not. In this
case, the fit should be rejected. An example of an acceptable and unacceptable description of Kd app is given
in Figure 4-7. The unacceptable fit on the right hand side of Figure 4-7 illustrates that the goodness of fit
cannot be assessed visually based on the mass and liquid phase concentrations alone.
Draft - April 2010
Page 18 of 40
Figure 4-7. Example of acceptable (left) and unacceptable (right) description of the apparent Kd values by the
aged sorption model
160
120
100
120
Mass (microg)
Mass (microg)
140
100
80
60
40
Model
20
80
60
40
20
Measurements
0
Model
0
0
20
40
60
80
100
0
20
Time (d)
1.0
0.5
0.0
20
40
60
Time (d)
80
Concentration (microg/l)
Concentration (microg/l)
Measurements
1.5
0
60
80
100
Time (d)
2.0
Model
40
Measurements
80
70
60
50
40
30
20
10
0
100
Model
0
20
40
60
Time (d)
Measurements
80
100
0.08
30
0.07
25
15
10
5
Model
Measurements
0
Kd app (L/kg)
Kd app (L/kg)
0.06
20
0.05
0.04
0.03
0.02
Model Kd
0.01
Measurements
0
0
20
40
60
Time (d)
80
100
0
20
40
60
Time (d)
80
100
The apparent Kd is not included in the fitting, a statistical goodness of fit criterion can thus not be calculated
and the assessment must be based on the visual agreement alone. However, the visual assessment of the
Kd app values is underpinned by the analysis of the fitted parameters. For example, the aged sorption
parameters fNE and kd are likely to be close to zero with large confidence intervals where there is no or only a
very small increase in Kd app.
In addition to the apparent Kd values, the simulated mass sorbed in the equilibrium and non-equilibrium
domain should be plotted against time. A robust fit is more likely when the mass in the non-equilibrium
domain shows a phase of decline during the experimental period. Fits are also more robust when nonequilibrium sorption is an important component of the whole system (i.e. the mass in the non-equilibrium
domain must not be negligible compared to the mass in the equilibrium domain). Examples are given in
Appendix 1.
Please note that using the best combination of parameters does not guarantee a good fit to the measured
data. If the model is not appropriate to describe measured behaviour, even the best possible parameter
combination for that model will not give an adequate fit to the data. The model will not be able to describe
the data, for example, if degradation is biphasic for reasons other than time-dependent sorption, or if
degradation shows a lag-phase. Always evaluate the visual fit to decide if a model is acceptable.
Draft - April 2010
Page 19 of 40
Chi2-test
4.5.2
FOCUS (2006) proposed a χ -test to evaluate the goodness of fit of degradation kinetics. A modified version
of the test should be applied to aged sorption data:
2
t
(Pi − Oi )2
i =1
( err / 100 x Oi )2
χ2 = ∑
(16)
The calculated χ for a specific fit may be compared to tabulated χ f ,α
2
2
where
t
= number of time points for mass plus number of time points for concentration
Pi = predicted value for measurement i
Oi = observed value for measurement i (replicates must be averaged to give a single value for each time
point)
err = measurement error percentage
f
= degrees of freedom = t minus number of model parameters
α = probability that one may obtain the given or higher χ2 by chance.
The χ -test considers the deviations between observed and predicted values relative to the uncertainty of the
measurements. Ideally, the measurement uncertainty at each time point should be determined from
numerous replicate values. Such replicate values are rarely available. Therefore, a pragmatic approach to
define the measurement variation was proposed by FOCUS (2006). The error of the measurements was
simply defined as a percentage of the average of all measurements. This implies that the absolute error is
identical for all measurements (i.e. for all time points). This is consistent with the recommendation of
unweighted fitting by FOCUS (2006). In contrast, the guidance on aged sorption studies proposes fitting to
weighted data for mass and liquid phase concentrations. The definition of the error has been changed to
reflect the assumptions that underlie weighted fitting. The error is now defined as a percentage of each
individual measurement (see denominator in Equation 16). As a result, the relative error is the same for all
measurements (i.e. all concentrations can be measured with the same relative precision). The absolute error
is now larger for large measurements.
2
The χ -test can be used to test the agreement between calculated and observed for a given fit. A suitable
model should pass the test at a significance level of 5%. However, this assessment is only possible if the
percent error is known. This is often not the case. Instead, the minimum error percentage at which the test is
2
passed (i.e. where the calculated value of χ is equal to or smaller than the standard tabulated value at the
5% significance level and the given degrees of freedom) can be directly derived from Equation 17.
2
err =
1
χ 2 tabulated
m
(Pi − Oi )2
i =1
Oi
∑
(17)
2
χ2tabulated = standard tabulated value at the 5% significance level and the given degrees of freedom
FOCUS (2006) recommends to calculate a χ error value for parent compounds and for metabolites
separately although the data for both compounds are fitted in a single optimisation. This division is
necessary because unweighted fitting is carried out and because the parent and metabolite data differ in
2
magnitude. The modified definition of the error in the χ test for aged sorption studies allows to calculate a
2
single χ error value for the mass and aqueous phase concentrations.
2
It is currently not possible to set a value for the error percentage as a criterion for acceptable fits to aged
sorption data. The value of 15% that is often used as a guide for degradation kinetics must not be applied to
aged sorption studies. More experience with fitting of the two-site sorption model is required before
recommendations can be made. Until specific guidance on the acceptable error value becomes available,
only the visual fit and relative standard error (see below) should be used as criteria for goodness of fit. It is,
2
however, useful to compare the χ error for the first-order equilibrium model (Section 4.4.2) and the aged
sorption model. The use of an aged sorption model is only justified if the error percentage is smaller. Again,
a minimum difference cannot be specified at this stage.
Draft - April 2010
Page 20 of 40
4.5.3
Confidence intervals and relative standard error
A confidence interval is an estimate of the uncertainty in a model parameter. The underlying assumption is:
If the experiment and the estimation procedure are repeated infinitely often, then the true value of the
parameter lies within the confidence interval with the chosen probability. The narrower the confidence
interval, the greater the precision with which the parameter can be estimated. Wide confidence intervals can
be caused by correlation between parameter values, parameter insensitivity, variability in the data, or the fact
that the model cannot describe the data.
TM
Optimisation tools such as ModelMaker
or PEST (used for optimising PearlNeq) give the optimised
parameters values together with the standard error or 95% confidence interval for each optimised parameter.
The standard error and the confidence interval should be converted into a relative standard error (RSE) as
follows:
RSE =
Standard error
υ
RSE =
95% Confidence Interval
4υ
(16,17)
where υ is the fitted parameter value. The confidence interval (upper limit minus lower limit) is divided by a
factor 4 to calculate the estimated standard deviation (or standard error) of the parameter fit. This is because
the width of the 95% confidence interval equals 4 times the standard deviation based on a normal
distribution (the fitted value plus or minus 2 × the standard deviation).
Wide confidence intervals imply that the parameters are very uncertain. Where 0 is included in the
confidence interval, there is not enough evidence that non-equilibrium sorption is a significant process. It is
difficult to set clear cut-off criteria for acceptable confidence intervals and relative standard errors. Based on
an analysis by DEFRA (2010), it is proposed that the RSE for any of the fitted parameters should not be
greater than 0.25. This implies that the width of the 95% confidence interval must not be greater than 100%
(i.e. ±50% of the parameter estimate).
4.5.4
Additional acceptance criteria
In addition to the criteria above, the following conditions must be met:
•
•
•
•
5
The fitted KOM,EQ value must be within ± 20% of the batch value from the same soil.
The fitted initial mass must be within ± 15% of the measured initial amount. If the mass at the time of
treatment is not measured, then the fitted initial mass can be compared with the added amount. It
should, however, be noted that the discrepancy is then also dependent on the experimental recovery
inherent in the extraction method.
The parameter fNE must be >0.001 and <10.
The parameter kd must be >0.00001 and <0.5.
Use of aged sorption parameters in regulatory exposure assessments
A sensitivity analysis by DEFRA (2010) showed that pesticide leaching models can be very sensitive for
changes in aged sorption parameters. It is thus very important to use robust parameter values in regulatory
exposure assessments.
5.1
Sources of input data for regulatory exposure assessments
5.1.1
Estimation of sorption parameters from soil or pesticide properties
Research by van DEFRA (2009) and DEFRA (2010) demonstrated that parameters of a two-site aged
sorption model can be very variable. Aged sorption parameters for the same pesticide can differ strongly
between different soils. There is no clear relationship between aged sorption parameters and soil or
pesticide properties. This was confirmed by an analysis by Sur et al. (2009). This is partly due to the
correlation between the two key parameters of the model fNE and kd. Different combinations of the two
parameters can give a similar result. A statistical relationship between a single parameter and soil or
pesticide properties is thus difficult to establish.
It is thus not recommended to estimate fNE and kd for new pesticides from soil or pesticide properties.
Draft - April 2010
Page 21 of 40
5.1.2
Default values
The large variation of the aged sorption parameters between studies implies that it is very difficult to
recommend default values for fNE and kd for use in the modelling of pesticide leaching. The increase in
sorption over time can be small for some pesticide-soil combinations. The use of default values would underestimate pesticide concentrations in leachate in these cases. Aged sorption parameters must instead be
determined in experiments with a range of soils.
5.1.3
Experimental laboratory incubation studies
Laboratory incubation studies are the recommended method to derive aged sorption parameters. Aged
sorption should be measured in four soils, ideally as part of OECD 307 laboratory degradation studies. Batch
sorption data are also needed for each of the four soils. Where possible, the soils used in the
degradation/aged sorption study should be identical with those used in the standard regulatory batch
sorption studies. Additional requirements for the methodology of aged sorption incubation studies are
outlined in Section 3.
DEFRA (2010) showed that the fitting of the two-site aged sorption model is less robust for pesticides with
small fNE and kd values. This is not expected to cause many problems in regulatory exposure assessments
because time-dependent sorption has only a very small effect on simulated concentrations in leachate for
these compounds. The implemetation of aged sorption as a higher tier option is thus not very meaningful for
pesticides with small fNE and kd values and robust parameters are not required. The model fitting is also less
robust for substances with intermediate fNE and kd values and long DegT50 values or for weakly sorbed
compounds. These compounds can be the cause of concern in the lower tier leaching assessments and it
may be desirable to consider aged sorption at the higher tiers. It is thus important to take measures that
reduce the likelihood of rejecting a fit of the aged sorption model for compounds with long half-lives and
small sorption coefficients. For example, it may be beneficial to include additional time points or replicates.
Prolonging the experimental period will also be beneficial, within the 120-day limit.
It must be noted that the optimisation of aged sorption parameters from laboratory aged sorption studies
does not replace the estimation of regulatory trigger values and degradation endpoints for modelling
according to the guidance given by FOCUS (2006). DegT50 values are derived in the aged sorption study
whereas DT50 values are required as regulatory persistence endpoints. Although first-order degradation
endpoints could be estimated from DegT50 values, this is not recommended. This is because the guidance
given by FOCUS (2006) differs from the guidance given in this document (e.g. comparison of first-order fit
with other kinetics, unweighted fitting, different goodness of fit criteria).
5.1.4
Column and field studies
As outlined by DEFRA (2010), it is not recommended to determine aged sorption parameters for use in
regulatory leaching modelling from column or field studies. Column and field studies can, however, be
presented as additional evidence that aged sorption is a relevant process for the substance of interest. For
example, it could be demonstrated that:
• The availability for leaching from columns declines over time. The decline is stronger than expected
from degradation alone.
• The apparent Kd value calculated from the field data increases over time
5.2
Aged sorption in the tiered pesticide leaching assessment
Regulatory leaching assessments are conducted in several tiers. At the lower tier, sorption is assumed to be
at equilibrium and degradation follows first-order kinetics (FOCUS, 2000). Sorption parameters that are used
as input for regulatory leaching models are taken from standard batch adsorption studies with at least five
soils. The average Koc or Kom value and the average Freundlich exponent (median if the number of soils is
large) are entered into the leaching model. Degradation endpoints for modelling are determined according to
FOCUS (2006) from standard aerobic soil degradation studies. The geometric mean of degradation
endpoints in at least four soils is used for modelling.
Time-dependent sorption can be considered within the regulatory procedure as a higher tier option. Aged
sorption parameters for parent compounds derived with PEARLNEQ can be used directly in the pesticide
leaching model FOCUS PEARL (all versions), PELMO (FOCUS GWII version) and PRZM (FOCUS GWII
version). For use in MACRO 5.0 onwards, the parameters must be converted (Equation 9 and 11).
Draft - April 2010
Page 22 of 40
In this section, it is discussed how the data from the aged sorption studies could be used at the higher tier of
the regulatory assessment and how lower and higher tiers of the assessment could be integrated. The
issues that need to be considered are:
•
It is recommended to carry out a minimum of four aged sorption studies. The question arises
whether the parameters from the four studies should be averaged or used individually. If individual
model runs are undertaken for each parameter combination, then guidance must be given on
th
whether the resulting 80 percentile concentrations in leachate should be averaged or not. The
derived parameters may not be reliable for all four experiments. It is not clear at this stage how this
problem should be addressed in the regulatory leaching assessment.
•
It must be considered whether the model input parameters on degradation and equilibrium sorption
should only be taken from the higher tier aged sorption studies or include the endpoints used at the
lower tier. If only the higher tier data are used, then information from the lower tier laboratory studies
and information from field dissipation studies would be overlooked. It is proposed to carry out four
aged sorption studies. The soils selected for these experiments may differ from those used in the
lower tier experiments. Degradation and equilibrium sorption parameters derived from the aged
sorption studies could, by chance, differ considerably from the endpoints used at the lower tier. The
higher tier leaching assessment could, thus, indicate a greater potential risk for leaching than the
lower tier even though aged sorption is included. It is, therefore, important to identify the conditions
that need to be met before the higher tier assessment overrides the lower tier.
An alternative option is to average the input parameters on degradation and equilibrium sorption
3
over the higher and lower tier studies . However, aged sorption parameters are highly correlated
with the parameters on degradation and equilibrium sorption. Each study gives a unique set of
parameters. Combining fNE and kd values from the aged sorption studies with average degradation
and equilibrium sorption parameters that include data from different soils is thus problematic.
The proposed procedure is outlined below:
It is not recommended to average aged sorption parameters prior to the leaching modelling. It is also not
appropriate to take the parameters fNE and kd from the aged sorption study and combine these with
degradation or batch sorption equilibrium values from independent studies. This is because fNE and kd are
not only closely correlated with each other, but also with the other parameters of the aged sorption model.
Several runs with the leaching model should thus be undertaken for each FOCUS leaching scenario: one run
for each combination of parameters derived from the incubation studies (fNE, kd, KOM,EQ, DegT50) and the
Freundlich exponent used in the optimisation. DegT50 values must be standardised to a moisture of pF 2
prior to modelling based on the guidance by FOCUS (2000 and 2002), unless the study was undertaken at
this moisture. It is assumed that the correction factors for DegT50 values are the same as those for firstorder DT50 values.
The proposed guidance implies that degradation and equilibrium sorption endpoints that are used for the
lower tier leaching assessment are not considered in the higher tier assessment. The higher tier leaching
simulations are based solely on the parameters derived from the aged sorption studies (unless the
parameters are not acceptable, see below). It is suggested that the simulated potential for movement to
groundwater at the higher tier overrides the lower tier result. This may not be justified where degradation and
equilibrium sorption in the aged sorption study are by chance very different from those in the lower tier
studies. This can arise because both types of experiments are carried out with only a small number of soils.
The decision on whether the higher tier or lower tier takes precedence should be based on a comparison of
4
the endpoints. If there is no systematic difference between degradation and equilibrium sorption parameters
from the two sets of studies, then the leaching assessment based on time-dependent sorption should
override the lower tier leaching assessment. The criteria for an acceptable difference between the endpoints
will have to be specified before this approach can be applied in practice.
th
th
The 80 percentile concentration in leachate at 1-m depth should be recorded for each simulation. An 80
percentile concentration must be calculated for each combination of parameters derived from the incubation
3
Note that the degradation endpoints from a lower tier degradation study must be converted to DegT50
values prior to use in simulations that account for aged sorption, see FOCUS (2010).
4
Note that the parameter DegT50 cannot be used in this comparison. The first-order DT50 must instead be
fitted against the total mass measured in the aged sorption study (this is already part of the guidance for
aged sorption studies, so no additional effort is required).
Draft - April 2010
Page 23 of 40
studies (fNE, kd, KOM,EQ, DegT50). Sorption is set to equilibrium sorption only for soils where the fitting of the
aged sorption model was not successful. The results for all soils are then averaged.
•
If the model fitting leads to an acceptable parameter combination, based on the criteria outlined in
Section 4, then these parameters should be used in the modelling.
•
If the model parameters are rejected for one or more soils, then the 80 percentile concentrations in
th
leachate for all of these cases are substituted by the 80 percentile concentration from the lower
tiers of the regulatory assessment. This will be a single value for each of the nine FOCUS leaching
scenarios, calculated assuming first-order degradation and equilibrium sorption. It may be based on
degradation and sorption endpoints that were derived from studies which were carried out prior to
the aged sorption experiments.
th
The procedure is illustrated for an example below:
Aged sorption studies were undertaken with four soils. Acceptable aged sorption parameters are available
th
for soils A-C. These are entered into FOCUS PEARL and 80 percentile concentrations in leachate are
recorded (Table 5-1) for e.g. the FOCUS Hamburg scenario. The parameters for soil D failed the criteria and
the lower tier value is used instead in the average calculation.
Table 5-1. Example results of a higher tier leaching assessment accounting for aged sorption
Soil
Parameters
acceptable?
Yes
Yes
Yes
No
A
B
C
D
Leaching modelling based on
Aged sorption
Aged sorption
Aged sorption
First-order degradation and equilibrium
sorption (lower tier endpoints averaged
over several soils)
Average
5.3
th
-1
80 percentile concentration (µg L ) for
FOCUS Hamburg scenario
0.050
0.089
0.031
0.254
0.106
Special considerations for metabolites
The fate of metabolites can be investigated in soils treated with the parent compound. Alternatively, they can
be added directly to the soil. It is not recommended to estimate aged sorption parameters for metabolites
that are formed during degradation of the parent. Fitting the aged sorption model simultaneously to the data
for the parent compound and metabolite involves optimising ten parameters (Mp,ini, fNE, kd, KOM,EQ, DegT50 for
the parent compound and fNE, kd, KOM,EQ, DegT50 and the formation fraction of the metabolite). This is very
challenging and it is likely that the aged sorption parameters for the metabolite fail the acceptability criteria
outlined in this document.
Aged sorption parameters from studies with the metabolite as the added substance can, in principle, be
derived based on the proposed guidance. There is, however, no evidence that the sorption behaviour of the
metabolite is the same when it is formed slowly and continuously from the parent or applied at once to the
soil.
There are a number of additional issues that need to be considered when aged sorption parameters for
metabolites are used in regulatory leaching assessments:
•
•
If the leaching of the parent and metabolite is calculated simultaneously, a formation fraction for the
metabolite must be entered into the model. This cannot be derived from the aged sorption study and
must be obtained from a degradation study with the parent compound as the added substance. It is
not clear whether it is valid to combine aged sorption parameters from a study with the metabolite as
the added substance with a formation fraction from an experiment with the parent as the applied
compound.
If the parent compound is also subject to aged sorption, then it must be ensured that the parameters
for the parent and metabolite are consistent (e.g. it may be necessary to request that these are
generated with the same soil).
Draft - April 2010
Page 24 of 40
6
References
Boesten, J.J.T.I., Tiktak, A. and van Leerdam R.C. (2007). Manual of PEARLNEQ v4. ALTERRA,
Wageningen: WOT Natuur & Milieu (Workdocuments 71) 34 pp.
Cox, L. and Walker, A. (1999). Studies of time-dependent sorption of linuron and isoproturon in soils.
Chemosphere 38:2707-2718.
DEFRA (2004). Time-dependent sorption processes in soil. Report for DEFRA project PS2206. Warwick
HRI. http://randd.defra.gov.uk/Document.aspx?Document=PS2206_3831_FRP.doc
DEFRA (2009). Characterisation and modelling of time-dependent sorption of pesticides. Report for DEFRA
project PS2228. The Food and Environment Research Agency.
http://randd.defra.gov.uk/Document.aspx?Document=PS2228_7878_FRP.doc
DEFRA (2010). Development of guidance on the implementation of aged soil sorption studies into regulatory
exposure assessments. Report for DEFRA project PS2235.The Food and Environment Research
Agency.
th
Doherty (2005). PEST. Model-independent parameter estimation. 5 edition. Watermark Numerical
Computing. www.sspa.com/pest, version 9.01. http://www.pesthomepage.org/files/pestman.pdf
FOCUS (2000). FOCUS groundwater scenarios in the EU pesticide registration process. Report of the
FOCUS Groundwater Scenarios Workgroup, EC Document Reference Sanco/321/2000 rev 2. 202pp.
FOCUS (2002). Generic Guidance for FOCUS groundwater scenarios, Version 1.1. 61 pp.
FOCUS (2006). Guidance document on estimating persistence and degradation kinetics from environmental
fate studies on pesticides in EU registration, report of the FOCUS work group on degradation kinetics,
EC document reference Sanco/10058/2005 Version 2.0, 434 pp.
FOCUS (2010). Assessing Potential for Movement of Active Substances and their Metabolites to Ground
Water in the EU. Report of the FOCUS Ground Water Work Group, EC Document Reference
Sanco/???/2010 Version 1, 594 pp.
Gurney, A.J.R. and Hayes, S.E. 2007. Non-equilibrium sorption and degradation of pesticides in soil:
analysis of laboratory aged sorption data using ModelMaker. In: A.A.M. Del Re, E. Capri, G.
Gragoulis, and M. Trevisan, eds. Environmental Fate and Ecological Effects of Pesticides
(Proceedings of the XIII Symposium Pesticide Chemistry). La Goliardica Pavese, Pavia Italy, ISBN
978-88-7830-473-4, pp. 245-253.
Leistra, M., Van der Linden, A.M.A., Boesten, J.J.T.I., Tiktak, A. and Van den Berg, F. (2001). PEARL model
for pesticide behaviour and emissions in soil-plant systems. Description of processes. Alterra report
013, Alterra, Wageningen, RIVM report 711401009, Bilthoven, The Netherlands. Available at
http://www.pearl.pesticidemodels.eu.
Larsbo, M. and Jarvis, N. (2003). MACRO 5.0. A model of water flow and solute transport in macroporous
soil. Technical description. Report. Department of Soil Science. Swedish University of Agricultural
Sciences. 49 pp.
MatLab (2007). MatLab Version 7.4.0.287 (R2007a), Optimisation Toolbox, Statistics Toolbox, MatLab
Compiler. The MathWorks Inc., USA. www.mathworks.com
OECD (2000). OECD Guidelines for the Testing of Chemicals. Test No 106: “Adsorption-Desorption using a
Batch Equilibrium Method” - Organization for Economic Cooperation and Development.
OECD (2002). OECD Guidelines for the Testing of Chemicals. Test No 307: Aerobic and anaerobic
transformation in soil. Organization for Economic Cooperation and Development.
Streck, T., Poletika N.N., Jury, W.A. and Farmer,W.J. (1995). Description of simazine transport with ratelimited,two-stage, linear and nonlinear sorption. Water Resources Research 31:811-822.
Sur, R., Menke, U., Dalkmann, P., Paetzold, S., Keppler, J. and Goerlitz, G. (2009). Comparative evaluation
of time-dependent sorption data of pesticides. Proceedings of the Conference on Pesticides in Soil,
Water and Air, 12-14 September 2009, York, UK.
Walker, A. and Jurado-Exposito, M. (1998). Adsorption of isoproturon, diuron and metsulfuron-methyl in two
soils at high soil:solution ratios. Weed Research 38:229-238.
Wauchope, R.D., Yeh, S., Linders, J.B.H.J., Kloskowski, R., Tanaka, K., Rubin, B., Katayama, A, Kordel, W.,
Gerstl, Z., Lane, M. and Unsworth J.B. (2002). Pesticide soil sorption parameters: theory, measurement,
uses, limitations and reliability. Pest Management Science 58:419-445.
Draft - April 2010
Page 25 of 40
Appendix 1: Fitting of a two-site aged sorption model with PEARLNEQ to two
example datasets
Example 1
An aged sorption laboratory incubation study was carried out using the experimental design described in
Section 3. The experimental conditions are shown in Table A1-1 and the measurements are given in Table
A1-2.
Table A1-1. Experimental conditions of the laboratory aged sorption study (example 1)
Parameter
Applied mass of pesticide
Mass of dry soil
Moisture
Water added for desorption
OC
OM
KF,om, batch
Freundlich exponent batch
Temperature
Limit of quantification in soil
Limit of quantification in CaCl2
Unit
µg
g
mL
mL
%
%
-1
mL g
o
C
µg
-1
µg mL
Value
20
8.52
1.48
20
1.47
2.53
246
0.830
20
4.00
0.026
Table A1-2. Measured data and calculated sorption and apparent Kd values (example 1)
Sampling
time
(days)
0
0
0
1
1
1
3
3
3
7
7
7
14
14
14
28
28
28
43
43
43
57
Total mass in organic
solvent extract
(µg)
20.179
20.401
20.088
20.293
20.307
20.381
19.192
19.123
18.931
18.737
18.575
18.235
17.490
17.596
17.853
16.233
16.202
16.263
14.929
14.993
15.231
13.846
Concentration in
extraction solution
-1
(µg mL )
0.2346
0.2304
0.2321
0.2243
0.2231
0.2212
0.1830
0.1871
0.2009
0.1843
0.1831
0.1780
0.1678
0.1647
0.1632
0.1295
0.1287
0.1271
0.1128
0.1083
0.1089
0.0947
57
57
71
71
71
82
82
82
13.779
13.712
13.613
13.165
12.796
12.489
12.423
11.930
0.0911
0.0966
0.0850
0.0821
0.0896
0.0799
0.0792
0.0793
Draft - April 2010
Sorbed amount
Apparent Kd
µg
mL g
1.78
1.81
1.77
1.82
1.82
1.83
1.79
1.77
1.72
1.73
1.72
1.69
1.63
1.65
1.68
1.58
1.58
1.59
1.47
1.49
1.51
1.39
7.57
7.87
7.64
8.10
8.16
8.29
9.79
9.48
8.54
9.41
9.38
9.50
9.71
10.02
10.32
12.19
12.25
12.50
13.01
13.73
13.90
14.64
1.39
1.37
1.38
1.34
1.28
1.26
1.26
1.20
15.24
14.14
16.28
16.30
14.24
15.81
15.88
15.14
-1
Page 26 of 40
Sorption was calculated from the measurements as:
Sorbed amount [ µg g −1 ] =
Total mass [ µg ] − (Liquid phase concentration [ µg mL−1 ] × Volume of liquid [mL])
Mass of soil [ g ]
The apparent Kd (Kd app) was calculated as:
K d app [ µg mL−1 ] =
Mass sorbed [ µg g −1 ]
Liquid phase concentration [ µg mL−1 ]
Measurements on day 0 and day 1 were discarded. All values are above the LOQ.
In order to carry out the non-equilibrium parameter estimation procedure in PEARLNEQ, the .mkn file of the
PEARLNEQ package has to be compiled following the instructions in the PEARLNEQ manual. The make file
of PEARLNEQ for the example case is shown below. Note that the current version of PEARLNEQ 4.0 does
not allow the fitting of replicates, a modified research version was used for this purpose.
The starting value for the initial mass (19.55 µg) and DegT50 (117.61 days) were derived by fitting a firstorder model to the data. The starting value for KOM,EQ was set to the batch value. Four starting value
combinations were tested for fNE and kd.
PEARLNEQ .mkn file for example case 1
* Model control
Yes
0.0
120.0
0.01
ScreenOutput
TimStart
TimEnd
DelTim
(d)
(d)
(d)
Start time of experiment
End time of experiment
Time step of Euler's integration procedure
* System characterisation
19.55
MasIni
8.52
MasSol
1.48
VolLiqSol
20.0
VolLiqAdd
0.0253
CntOm
(ug)
(g)
(mL)
(mL)
(kg.kg-1)
Initial guess of initial mass
Mass of soil in incubation jar
Volume of liquid in the moist soil
Volume of liquid ADDED
Organic matter content
* Sorption parameters
1.0
ConLiqRef
0.830
ExpFre
246
KomEql
0.2
FacSorNeqEql
0.004
CofRatDes
(mg.L-1)
(-)
(L.kg-1)
(-)
(d-1)
Reference liquid concentration
Freundlich exponent
Coefficient for equilibrium sorption
Initial guess of ratio KfNeq/KfEql
Initial guess of desorption rate constant
* Transformation
117.61
20.00
110.00
(d)
(C)
(kJ.mol-1)
Initial guess of half-life at ref temperature
Reference temperature
Initial guess of molar activation energy
parameters
DT50Ref
TemRefTra
MolEntTra
* Temperatures at which the incubation experiments have been carried out
table Tem (C)
1 20
end_table
* Provide the results
* Tim
Tem
* (d)
(C)
table Observations
3
20
19.192
3
20
19.123
3
20
18.931
7
20
18.737
7
20
18.575
7
20
18.235
14
20
17.490
14
20
17.596
14
20
17.853
28
20
16.233
28
20
16.202
28
20
16.263
43
20
14.929
43
20
14.993
43
20
15.231
Draft - April 2010
of the measurements
Mass
(ug)
ConLiq
(ug/mL)
0.1830
0.1871
0.2009
0.1843
0.1831
0.1780
0.1678
0.1647
0.1632
0.1295
0.1287
0.1271
0.1128
0.1083
0.1089
Page 27 of 40
57
20
13.846 0.0947
57
20
13.779 0.0911
57
20
13.712 0.0966
71
20
13.613 0.0850
71
20
13.165 0.0821
71
20
12.796 0.0896
82
20
12.489 0.0799
82
20
12.423 0.0792
82
20
11.930 0.0793
end_table
* Option for weights of observations:
*'equal' gives equal weights to all measurements
*'inverse' gives weight equal to inverse value of each measurement (if measurement is zero then
weight is 1.0)
inverse Opt_weights
Running the PEARLMK program produces a series of files that are necessary to run the PEST
optimisation. The key file is the PEST control file with the extension “*.PST”. The control file of the
example is shown below:
PEST control file for example case 1
pcf
* control data
restart
5 48
5 0
1 3 single point
5.0 2.0 0.1 0.01 15
3.0 4.0 1.0e-3
0.1
50 0.001 5 10 0.001 4
1 1 1
* group definitions and derivative data
FSNE
relative 0.01 0.00001 always_3
2.0 best_fit
CRD
relative 0.01 0.00001 always_3
2.0 best_fit
DT50
relative 0.01 0.00001 always_3
2.0 best_fit
MASINI relative 0.01 0.00001 always_3
2.0 best_fit
KOMEQ relative 0.01 0.00001 always_3
2.0 best_fit
* parameter data
FSNE
none factor
0.2
0.001
10.0 FSNE
CRD
none factor
0.004
0.00001
0.5 CRD
DT50
none factor
117.61
1.0
500.0 DT50
MASINI none factor
19.55
0.1
1000.0 MASINI
KOMEQ none factor
246
0.1
40000.0 KOMEQ
* observation data and weights
O1
19.192 0.05210506
O2
18.737 0.053369176
O3
17.490 0.057176527
O4
16.233 0.061602817
O5
14.929 0.066984902
O6
13.846 0.072221211
O7
13.613 0.07345857
O8
12.489 0.080070795
O9
0.1830 5.465808563
O10
0.1843 5.426927118
O11
0.1678 5.959860639
O12
0.1295 7.720238523
O13
0.1128 8.861637214
O14
0.0947 10.56126183
O15
0.0850 11.76580983
O16
0.0799 12.50794457
O17
19.123 0.052294311
O18
18.575 0.053835594
O19
17.596 0.056830513
O20
16.202 0.061720352
O21
14.993 0.066696807
O22
13.779 0.072572737
O23
13.165 0.075956964
O24
12.423 0.080494218
O25
0.1871 5.344996603
O26
0.1831 5.460953069
O27
0.1647 6.071034367
O28
0.1287 7.769100473
O29
0.1083 9.233575211
O30
0.0911 10.9825065
O31
0.0821 12.17780069
O32
0.0792 12.61950108
Draft - April 2010
1.00
1.00
1.00
1.00
1.00
0.00
0.00
0.00
0.00
0.00
Page 28 of 40
O33
18.931 0.052824376
O34
18.235 0.054840998
O35
17.853 0.056012089
O36
16.263 0.061489076
O37
15.231 0.065656019
O38
13.712 0.072929498
O39
12.796 0.07814921
O40
11.930 0.083820865
O41
0.2009 4.978075535
O42
0.1780 5.618356373
O43
0.1632 6.126011669
O44
0.1271 7.87036175
O45
0.1089 9.184578133
O46
0.0966 10.35478581
O47
0.0896 11.15958672
O48
0.0793 12.61596801
* model command line
..\..\neq_bin\PearlNeq example
* model input/output
example.tpl
example.neq
example1.ins
example.out
example2.ins
example.out
example3.ins
example.out
After PEST is started, PEST runs the PEARLNEQ model. This produces an output file as shown below. The
results of the output file are then compared to the measured data by the PEST program and the parameters
are changed until the sum of squared residues is minimised or the termination criteria specified in the pest
control file are met.
Output file for example case 1 (starting value combination 1)
*
*
*
*
*
*
*
*
*
-----------------------------------------------------------------------------Results from PEARLNEQ (c) MNP/RIVM/Alterra
PEARLNEQ version 4.1
PEARLNEQ created on 04-Jan-2010
Run ID
: example
Input file generated on
: 04-03-2010
------------------------------------------------------------------------------
-------------------------------------------------------------* System properties
* Mass of dry soil (g)
:
8.5200
* Volume of water in moist soil (mL)
:
1.4800
* Volume of water added (mL)
:
20.0000
* Initial mass of pesticide (ug)
:
19.3909
* Reference concentration (ug.mL-1)
:
1.0000
* Equilibrium sorption coeff (mL.g-1)
:
6.5461
* Non-equili. sorption coeff (mL.g-1)
:
2.8125
* Freundlich exponent (-)
:
0.8300
* Desorption rate coefficient (d-1)
:
0.0237
* Half-life transformation (d)
:
98.4763
* Reference temperature (K)
: 293.1500
-------------------------------------------------------------*
*
Temp
(C)
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
Time
(d)
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
13.0
14.0
15.0
16.0
17.0
18.0
19.0
Draft - April 2010
Mas
(ug)
19.39086700
19.25551616
19.12241967
18.99150935
18.86271938
18.73598622
18.61124851
18.48844703
18.36752460
18.24842601
18.13109796
18.01548898
17.90154937
17.78923116
17.67848801
17.56927516
17.46154942
17.35526904
17.25039374
17.14688457
ConLiq
(ug.mL-1)
0.20523346
0.20127622
0.19744824
0.19374467
0.19016084
0.18669226
0.18333462
0.18008378
0.17693577
0.17388675
0.17093305
0.16807114
0.16529761
0.16260921
0.16000279
0.15747534
0.15502395
0.15264585
0.15033834
0.14809885
XNeq
(ug.g-1)
0.00000000
0.02227986
0.04366741
0.06419473
0.08389278
0.10279139
0.12091936
0.13830447
0.15497348
0.17095224
0.18626566
0.20093778
0.21499177
0.22845000
0.24133401
0.25366462
0.26546189
0.27674514
0.28753306
0.29784363
Page 29 of 40
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
21.0
22.0
23.0
24.0
25.0
26.0
27.0
28.0
29.0
30.0
31.0
32.0
33.0
34.0
35.0
36.0
37.0
38.0
39.0
40.0
41.0
42.0
43.0
44.0
45.0
46.0
47.0
48.0
49.0
50.0
51.0
52.0
53.0
54.0
55.0
56.0
57.0
58.0
59.0
60.0
61.0
62.0
63.0
64.0
65.0
66.0
67.0
68.0
69.0
70.0
71.0
72.0
73.0
74.0
75.0
76.0
77.0
78.0
79.0
80.0
81.0
82.0
83.0
84.0
85.0
86.0
87.0
88.0
89.0
90.0
91.0
92.0
93.0
94.0
95.0
96.0
97.0
98.0
99.0
Draft - April 2010
17.04470396
16.94381559
16.84418438
16.74577646
16.64855910
16.55250067
16.45757065
16.36373951
16.27097874
16.17926079
16.08855904
15.99884777
15.91010212
15.82229805
15.73541235
15.64942258
15.56430702
15.48004473
15.39661541
15.31399946
15.23217794
15.15113251
15.07084545
14.99129963
14.91247846
14.83436592
14.75694650
14.68020518
14.60412746
14.52869930
14.45390709
14.37973769
14.30617838
14.23321682
14.16084110
14.08903966
14.01780133
13.94711528
13.87697103
13.80735843
13.73826763
13.66968912
13.60161366
13.53403230
13.46693638
13.40031749
13.33416748
13.26847847
13.20324280
13.13845303
13.07410198
13.01018264
12.94668825
12.88361222
12.82094818
12.75868992
12.69683144
12.63536690
12.57429062
12.51359709
12.45328097
12.39333706
12.33376030
12.27454578
12.21568872
12.15718450
12.09902859
12.04121659
11.98374425
11.92660740
11.86980200
11.81332410
11.75716987
11.70133558
11.64581759
11.59061236
11.53571642
11.48112642
11.42683906
11.37285116
0.14592490
0.14381412
0.14176420
0.13977294
0.13783823
0.13595803
0.13413038
0.13235339
0.13062524
0.12894421
0.12730860
0.12571680
0.12416727
0.12265850
0.12118907
0.11975758
0.11836271
0.11700317
0.11567774
0.11438523
0.11312450
0.11189444
0.11069401
0.10952218
0.10837798
0.10726046
0.10616871
0.10510186
0.10405907
0.10303953
0.10204246
0.10106710
0.10011274
0.09917867
0.09826423
0.09736876
0.09649165
0.09563230
0.09479013
0.09396458
0.09315511
0.09236122
0.09158239
0.09081816
0.09006806
0.08933164
0.08860847
0.08789815
0.08720027
0.08651445
0.08584031
0.08517751
0.08452569
0.08388452
0.08325369
0.08263289
0.08202181
0.08142017
0.08082770
0.08024412
0.07966919
0.07910264
0.07854425
0.07799378
0.07745101
0.07691571
0.07638770
0.07586676
0.07535269
0.07484533
0.07434448
0.07384896
0.07336063
0.07287831
0.07240184
0.07193108
0.07146588
0.07100609
0.07055158
0.07010223
0.30769421
0.31710154
0.32608176
0.33465045
0.34282263
0.35061277
0.35803485
0.36510234
0.37182823
0.37822505
0.38430488
0.39007939
0.39555980
0.40075696
0.40568133
0.41034297
0.41475162
0.41891664
0.42284708
0.42655167
0.43003879
0.43331658
0.43639283
0.43927510
0.44197065
0.44448649
0.44682939
0.44900587
0.45102220
0.45288446
0.45459848
0.45616991
0.45760418
0.45890653
0.46008200
0.46113548
0.46207167
0.46289507
0.46361008
0.46422089
0.46473157
0.46514602
0.46546803
0.46570122
0.46584911
0.46591508
0.46590239
0.46581417
0.46565346
0.46542318
0.46512614
0.46476506
0.46434255
0.46386113
0.46332324
0.46273121
0.46208730
0.46139370
0.46065249
0.45986570
0.45903527
0.45816309
0.45725096
0.45630064
0.45531378
0.45429203
0.45323693
0.45215000
0.45103267
0.44988635
0.44871238
0.44751205
0.44628662
0.44503729
0.44376521
0.44247150
0.44115724
0.43982346
0.43847115
0.43710129
Page 30 of 40
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
20.0
100.0
101.0
102.0
103.0
104.0
105.0
106.0
107.0
108.0
109.0
110.0
111.0
112.0
113.0
114.0
115.0
116.0
117.0
118.0
119.0
120.0
Draft - April 2010
11.31915958
11.26576128
11.21265329
11.15983270
11.10729669
11.05504248
11.00306738
10.95136875
10.89994400
10.84879062
10.79790613
10.74728813
10.69693427
10.64684222
10.59700975
10.54743462
10.49811469
10.44904783
10.40023196
10.35166506
10.30334512
0.06965790
0.06921847
0.06878382
0.06835385
0.06792843
0.06750746
0.06709084
0.06667846
0.06627024
0.06586608
0.06546589
0.06506958
0.06467706
0.06428826
0.06390310
0.06352150
0.06314338
0.06276868
0.06239733
0.06202925
0.06166894
0.43571479
0.43431254
0.43289541
0.43146421
0.43001975
0.42856278
0.42709405
0.42561427
0.42412412
0.42262425
0.42111530
0.41959788
0.41807257
0.41653994
0.41500053
0.41345486
0.41190344
0.41034675
0.40878525
0.40721939
0.40564960
Page 31 of 40
The results of the optimisation are written into a file with the extension .rec. Running the PEST optimisation
for the example case yields the results below.
Results for example 1, starting value combination 1 (fNE = 0.2, kd = 0.004)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
0.429644
2.374104E-02
98.4763
19.3909
258.738
95% percent confidence limits
lower limit
upper limit
0.371338
0.487951
1.749243E-02
2.998964E-02
92.7527
104.200
19.0690
19.7128
249.546
267.930
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
2.46139E-02
Results for example 1, starting value combination 2 (fNE = 0.2, kd = 0.05)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
0.429636
2.374208E-02
98.4757
19.3909
258.738
95% percent confidence limits
lower limit
upper limit
0.371291
0.487981
1.749984E-02
2.998432E-02
92.7490
104.202
19.0690
19.7128
249.546
267.931
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
2.46138E-02
Results for example 1, starting value combination 3 (fNE = 1.5, kd = 0.004)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
0.429642
2.374070E-02
98.4765
19.3909
258.739
95% percent confidence limits
lower limit
upper limit
0.371311
0.487973
1.749748E-02
2.998392E-02
92.7500
104.203
19.0690
19.7128
249.547
267.931
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
2.46139E-02
Results for example 1, starting value combination 4 (fNE = 1.5, kd = 0.05)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
0.429632
2.374298E-02
98.4752
19.3909
258.738
95% percent confidence limits
lower limit
upper limit
0.371299
0.487965
1.749982E-02
2.998614E-02
92.7486
104.202
19.0690
19.7128
249.545
267.930
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
2.46138E-02
The four starting value combinations gave almost identical objective functions (sum of squared weighted
residuals = phi) and parameter values. Relative standard errors (RSE) were calculated from the confidence
intervals for starting combination 2 and 4 (those with the smallest phi) based on equation 16. Combination 2
gave marginally smaller RSE than combination 4 and this combination was chosen for further analysis.
Draft - April 2010
Page 32 of 40
Table A1-3. Optimisation results for example 1, starting combination 2
Parameter
fNE
kd
DegT50
MasIni
KOM,EQ
Optimised
value
0.4296
0.0237
98.48
19.39
258.8
RSE
0.068
0.131
0.029
0.008
0.018
RSE <0.25?
yes
yes
yes
yes
yes
Check
Result
Within ±20% of added mass (20 µg)?
-1
Within ±20% of batch value (246 mL g )?
yes
yes
The RSE values of all parameters are below 0.25. The fitted mass and KOM,EQ are within the acceptable
limits.
The visual fit of the aged sorption model is shown in Figure A1-1. The top graphs show the simulated mass
and concentrations in the liquid phase compared with the measured data. The residuals (deviations of
simulated minus observed) are presented in the middle graphs. The bottom left graph shows the apparent
Kd value compared with the values calculated from the measured data. The apparent Kd value was not
included in the model fitting and is only presented to support the interpretation of the results. The graph on
5
the bottom right shows the simulated sorbed mass in the equilibrium and non-equilibrium phase.
The visual fit to the mass and concentrations in the liquid phase is very good. The residuals are small and
randomly distributed around the zero line. The Kd app values show a clear increase in sorption over time and
are well described by the model. The sorbed mass in the non-equilibrium phase increases up to
approximately 60 days and starts to decline very slightly thereafter.
The aged sorption model was compared with an equilibrium model that does not account for time-dependent
sorption processes (Figure A1-2). This was achieved by fixing fNE and kd at zero. The model is not able to
describe the observed data. The residuals show systematic deviations from the zero line There is a large
discrepancy between the calculated Kd app values and the simulated line (the small increase in the simulated
Kd app values is due to non-linear sorption; Freundlich exponent = 0.83)
The χ test resulted in a very small error percentage (1.1%) for the aged sorption model. This is a
2
considerable improvement over the equilibrium model (χ error = 6.1%).
2
The overall conclusion on this example case is:
The fit of the aged sorption model to the data for example 1 is acceptable and the parameter values
can be used for modelling.
5
Note that the current version of PEARLNEQ (Version 4) provides output of the total mass, concentration in
liquid phase (ConLiq) and mass in the non-equilibrium phase (XNeq). The mass in the equilibrium phase
(Xeq) can be calculated as fitted KF,EQ x concentration in liquid phase ^ N.
Kd app = (Xeq+Xneq)/ConLiq.
Draft - April 2010
Page 33 of 40
Figure A1-1. Fitted vs measured mass and liquid phase concentrations and residuals for the aged sorption
model fitted to example 1
25
0.25
Measurements
20
Mass (microg)
Concentration (microg/l)
PEARLNEQ; Chi2=1.1
15
10
5
0
PEARLNEQ; Chi2=1.1
0.20
Measurements
0.15
0.10
0.05
0.00
0
20
40
60
80
100
0
20
Time (d)
Residuals conc.(microg/l)
1.5
1.0
0.5
0.0
-0.5
-1.0
-1.5
-2.0
80
100
80
100
0.04
0.03
0.02
0.01
0
-0.01
-0.02
-0.03
-0.04
0
20
40
60
80
100
0
20
40
Time (d)
60
Time (d)
2.0
18
16
14
12
10
8
6
4
2
0
Sorbed mass (microg/g)
Kd app (mL/g)
60
Time (d)
2.0
Residuals mass (microg)
40
PEARLNEQ
Calculated from measurements
Equilibrium sorption
1.5
Non-equilibrium sorption
1.0
0.5
0.0
0
20
40
Time (d)
Draft - April 2010
60
80
100
0
20
40
60
80
100
Time (d)
Page 34 of 40
Figure A1-2. Fitted vs measured mass and liquid phase concentrations and residuals for the equilibrium
sorption model fitted to example 1
25
Measurements
Mass (microg)
20
15
10
5
Concentration (microg/l)
0.25
Equilibrium model; Chi2=6.1
Equilibrium model; Chi2=6.1
0.20
Measurements
0.15
0.10
0.05
0.00
0
0
20
40
60
80
0
100
20
40
Time (d)
Residuals conc.(microg/l)
Residuals mass (microg)
2.0
1.5
1.0
0.5
0.0
-0.5
-1.0
-1.5
-2.0
80
100
80
100
0.04
0.03
0.02
0.01
0
-0.01
-0.02
-0.03
-0.04
0
20
40
60
80
100
0
20
40
Time (d)
60
Time (d)
2.5
18
16
14
12
10
8
6
4
2
0
Sorbed mass (microg/g)
Kd app (mL/g)
60
Time (d)
Equilibrium model
Calculated from measurements
Equilibrium sorption
2.0
Non-equilibrium sorption
1.5
1.0
0.5
0.0
0
20
40
Time (d)
Draft - April 2010
60
80
100
0
20
40
60
80
100
Time (d)
Page 35 of 40
Example 2
An aged sorption laboratory incubation study was carried out using the experimental design described in
Section 3. The experimental conditions are shown in Table A1-1 and the measurements are given in Table
A1-2.
Table A1-1. Experimental conditions of the laboratory aged sorption study (example 2)
Parameter
Applied mass of pesticide
Mass of dry soil
Moisture
Water added for desorption
OC
OM
KF,om, batch
Freundlich exponent batch
Temperature
Limit of quantification in soil
Limit of quantification in CaCl2
Unit
µg
g
mL
mL
%
%
-1
mL g
o
C
µg
-1
µg mL
Value
70
6.81
3.19
20
3.3
5.7
101
0.814
20
1.38
0.020
Table A1-2. Measured data and calculated sorption and apparent Kd values (example 2)
Sampling
time
(days)
0
0
0
1
1
1
3
3
3
7
7
7
14
14
14
28
28
28
42
42
42
56
Total mass in organic
solvent extract
(µg)
68.768
71.474
70.898
68.190
69.036
71.367
63.886
61.465
63.458
57.297
56.269
55.976
41.764
49.310
53.560
34.513
35.421
35.993
29.749
25.959
26.522
19.143
Concentration in
extraction solution
-1
(µg mL )
1.1157
1.1000
1.0949
1.0900
1.1122
1.1126
1.0157
0.9924
0.9906
0.8971
0.8654
0.8688
0.6786
0.6570
0.7042
0.4425
0.4679
0.4637
0.2861
0.2940
0.2986
0.2159
56
56
70
70
70
83
83
83
18.600
19.133
14.077
16.157
14.397
10.715
10.887
9.440
0.1926
0.1716
0.1313
0.1329
0.1132
0.0733
0.0786
0.0770
Sorbed amount
Apparent Kd
µg
mL g
6.30
6.75
6.68
6.30
6.35
6.69
5.92
5.65
5.94
5.36
5.32
5.26
3.82
5.00
5.47
3.56
3.61
3.71
3.39
2.81
2.88
2.08
5.65
6.14
6.10
5.78
5.71
6.01
5.83
5.69
6.00
5.97
6.14
6.06
5.63
7.62
7.76
8.05
7.71
7.99
11.86
9.56
9.64
9.61
2.08
2.23
1.62
1.92
1.73
1.32
1.33
1.12
10.78
12.96
12.33
14.44
15.26
18.07
16.93
14.60
-1
The aged sorption model was fitted to the mass and liquid phase concentration using PEARLNEQ. The
starting value for the initial mass (68.20 µg) and DegT50 (30.43 days) were derived by fitting a first-order
model to the data. The starting value for KOM,EQ was set to the batch value. Four starting value combinations
were tested for fNE and kd. The results are shown below:
Draft - April 2010
Page 36 of 40
Results for example 2, starting value combination 1 (fNE = 0.2, kd = 0.004)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
4.56555
3.924226E-04
27.0553
69.6762
108.445
95% percent confidence limits
lower limit
upper limit
-144.108
153.239
-1.309112E-02
1.387596E-02
25.7979
28.3127
66.8218
72.5306
98.9897
117.899
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
0.1938
Results for example 2, starting value combination 2 (fNE = 0.2, kd = 0.05)
Parameter
Estimated
value
fsne
10.0000
crd
1.770197E-04
dt50
27.0709
masini
69.6649
komeq
108.491
Objective function ----->
95% percent confidence limits
lower limit
upper limit
-663.370
683.370
-1.184228E-02
1.219632E-02
25.7824
28.3594
66.8081
72.5216
99.0147
117.967
Sum of squared weighted residuals (ie phi)
=
0.1935
Results for example 2, starting value combination 3 (fNE = 1.5, kd = 0.004)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
4.15434
4.330716E-04
27.0511
69.6616
108.359
95% percent confidence limits
lower limit
upper limit
-118.979
127.287
-1.299090E-02
1.385704E-02
25.7931
28.3091
66.8066
72.5167
98.9063
117.812
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
0.1939
Results for example 2, starting value combination 4 (fNE = 1.5, kd = 0.05)
Parameter
fsne
crd
dt50
masini
komeq
Estimated
value
10.0000
1.770196E-04
27.0709
69.6649
108.491
95% percent confidence limits
lower limit
upper limit
-660.714
680.714
-1.182025E-02
1.217429E-02
25.7823
28.3595
66.8082
72.5217
99.0156
117.967
Objective function ----->
Sum of squared weighted residuals (ie phi)
=
0.1935
Starting value combinations 2 and 4 resulted in fNE values that reached the upper limit allowed in the
optimisation (10). These fits are not acceptable. Starting value combinations 1 and 3 gave similar parameter
combinations. The sum of squared residuals was slightly smaller for combination 1 and this was chosen for
further analysis.
Table A1-5. Optimisation results for example 2, starting combination 1
Parameter
fNE
kd
DegT50
MasIni
KOM,EQ
Optimised
value
4.57
3.92E-04
27.06
69.68
108.45
RSE
RSE <0.25?
Check
Result
16.28
17.18
0.023238
0.020483
0.043592
no
no
yes
yes
yes
Within ±20% of added mass (70 µg)?
-1
Within ±20% of batch value (101 mL g )?
yes
yes
The RSE values of fNE and kd are well above 0.25 (Table A1-5). The fit is thus not acceptable. Figure A1-3
and Figure A1-4 illustrate the visual fit of the aged sorption model and the equilibrium model, respectively.
Note that graphs do not have to be generated when the fit is rejected, they are shown here to facilitate the
interpretation of the results.
Draft - April 2010
Page 37 of 40
The parameters fNE and kd are uncertain for two reasons:
•
The data on total mass are somewhat scattered.
•
The extent of non-equilibrium sorption that was observed within the experimental period is small for
this dataset.
There is a relatively small difference between the equilibrium model (Figure A1-4) and the aged
sorption model (Figure A1-3). The observed increase in Kd app values over time is only slightly larger
than can be explained by non-linear sorption on its own (as indicated by the discrepancy of the Kd app
line and symbols in Figure A1-4). The low importance of non-equilibrium sorption is also apparent
from the bottom right graph in Figure A1-3.
This example illustrates that a good visual agreement between measured and simulated data and a small χ
error value do not guarantee an acceptable fit. If the data are only weakly influenced by non-equilibrium
sorption, then the aged sorption model cannot be discriminated from the equilibrium model. The additional
parameters of the aged sorption model (fNE and kd) can then not be determined with sufficient confidence.
2
The overall conclusion on this example case is:
The fit of the aged sorption model to the data for example 1 is NOT acceptable and the parameter
values CANNOT be used for modelling.
Draft - April 2010
Page 38 of 40
Figure A1-3. Fitted vs measured mass and liquid phase concentrations and residuals for the aged sorption
model fitted to example 2
80
1.20
Mass (microg)
Concentration (microg/l)
PEARLNEQ; Chi2=3.6
70
Measurements
60
50
40
30
20
10
0
PEARLNEQ; Chi2=3.6
1.00
Measurements
0.80
0.60
0.40
0.20
0.00
0
20
40
60
80
100
0
20
Time (d)
Residuals conc.(microg/l)
8.0
4.0
0.0
-4.0
-8.0
-12.0
80
100
80
100
0.15
0.10
0.05
0.00
-0.05
-0.10
-0.15
0
20
40
60
80
100
0
20
40
Time (d)
60
Time (d)
8.0
20
18
16
14
12
10
8
6
4
2
0
Sorbed mass (microg/g)
Kd app (mL/g)
60
Time (d)
12.0
Residuals mass (microg)
40
PEARLNEQ
Calculated from measurements
7.0
Equilibrium sorption
6.0
Non-equilibrium sorption
5.0
4.0
3.0
2.0
1.0
0.0
0
20
40
Time (d)
Draft - April 2010
60
80
100
0
20
40
60
80
100
Time (d)
Page 39 of 40
Figure A1-4. Fitted vs measured mass and liquid phase concentrations and residuals for the aged sorption
model fitted to example 2
80
1.20
Mass (microg)
Concentration (microg/l)
PEARLNEQ; Chi2=7.4
70
Measurements
60
50
40
30
20
10
0
PEARLNEQ; Chi2=7.4
1.00
Measurements
0.80
0.60
0.40
0.20
0.00
0
20
40
60
80
100
0
20
Time (d)
Residuals conc.(microg/l)
8.0
4.0
0.0
-4.0
-8.0
-12.0
80
100
80
100
0.15
0.10
0.05
0.00
-0.05
-0.10
-0.15
0
20
40
60
80
100
0
20
40
Time (d)
60
Time (d)
8.0
20
18
16
14
12
10
8
6
4
2
0
Sorbed mass (microg/g)
Kd app (mL/g)
60
Time (d)
12.0
Residuals mass (microg)
40
PEARLNEQ
Calculated from measurements
7.0
Equilibrium sorption
6.0
Non-equilibrium sorption
5.0
4.0
3.0
2.0
1.0
0.0
0
20
40
Time (d)
Draft - April 2010
60
80
100
0
20
40
60
80
100
Time (d)
Page 40 of 40