A methodological approach to develop “contaminant migration–population effects” models

Ecological Modelling 220 (2009) 3280–3290

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

A methodological approach to develop “contaminant migration–population effects” models Luigi Monte ∗ ENEA CR Casaccia, via P. Anguillarese, 301, CP 2400, 00100, Rome, Italy

a r t i c l e

i n f o

Article history: Received 21 May 2009 Received in revised form 28 August 2009 Accepted 4 September 2009 Available online 30 September 2009 Keywords: Contaminant migration models Population dynamics Population dispersal Contamination of moving organisms Contamination and age-class effects

a b s t r a c t The main aim of the present work is to discuss the methodological approaches that underpin the “contaminant migration–population effects” models for the evaluation of the detriment to populations of moving organisms in environmental systems with spatial and time dependent pollution levels. A technique to couple the equations controlling the population dynamics and the pollutant dispersion is described and discussed. The domain of application and the limitations of the methodology are analysed and illustrated by some examples. Possible alternative approaches are briefly presented. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The assessment of the dispersion in the biosphere of toxic substances and of the consequent impacts on the biota populations is essential for the proper management of the aftermath of events of accidental introduction of pollutants into the environment. Recently this problem has more and more attracted the attention of radioecologists. In particular, radioecological studies (Beresford et al., 2008; Woodhead, 2003) are focusing on the development of methodological approaches and models for assessing the impact of high levels of radiation doses on the ecosystem in areas such as the Chernobyl Exclusion Zone and the Eastern Ural Radioactive Trace that were heavily contaminated by radioactive substances released following severe nuclear accidents (IAEA, 2006; Joint Norwegian-Russian Group, 1997; UNSCEAR, 2000; Sazykina and Kryshev, 2003a,b, 2006). It could be worthwhile to review and evaluate the modelling methodologies which can be used to deal with the above-mentioned problem not only for the management of the consequences of radioactive accidents but also for other environmental contamination emergencies. “Contaminant migration–population effects” models account for three groups of reciprocally interdependent processes:

∗ Tel.: +39 0630484645; fax: +39 0630483213. E-mail address: [email protected]. 0304-3800/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2009.09.001

(a) the dynamics and the dispersal of the biota populations; (b) the migration of the contaminants through the abiotic environmental components, the transfer from these to living organisms and, finally, the following dispersion of contaminant in the environment due to the biota movements; (c) the influence on the population dynamics of the detrimental effects caused by the environmental contamination. The web of mutual relations between these different processes is overwhelming complicated as biota populations, on the one hand, are affected by the environmental contamination and, on the other hand, contribute to the dispersion of contaminant in the biosphere (Goodwin et al., 2006; O’Toole et al., 2006). A great many models were developed to predict the transfer of pollutants through the environment (Håkanson and Peters, 1995; Håkanson and Bryhn, 2008; Jørgensen, 1978; Jørgensen, 1983; Jørgensen, 1984; Lauenroth et al., 1983; Scott, 2003; Smith and Comans, 1996; Smith et al., 2004) or to assess the movement of biota (Crist and Wiens, 1995; Li et al., 1996; Marsh and Jones, 1988; Turchin, 1991). Comparatively less attention has been paid to the comprehensive modelling of pollutant transport through the environment in connection with population dispersal and dynamics (Buonomo et al., 1999; Hallam et al., 1983). The main aim of the present work is to present and discuss the methodological approaches that underpin this last kind of models. We will describe: (a) some particular results that can be predicted by the models based on the methodology, such as the age-class effects on the distribution of contaminant in individual organisms of different age and the contamination dynamics in mov-

L. Monte / Ecological Modelling 220 (2009) 3280–3290

ing biota; and (b) the main limitation of the approach and, briefly, the possible alternative methods. The analysis of the foundations of the methodology is helpful to identify the application domains of the derived models. 2. Modelling the environmental dispersion of a contaminant transported by biota 2.1. Preliminary remarks We will start our analysis by reviewing the methods and discussing the hypotheses underpinning a general system of equations to predict the population dynamics and the transport of contaminants by living organisms through the environment. An equation controlling the concentration of a contaminant in moving biota can be derived by applying the so-called particle˜ tracking or “Lagrangian” method (Keats et al., 2007; Perianez, 2007; Rajar et al., 1997). Each individual is followed as it moves through space and the contaminant concentration in the organism is determined by accounting for the mechanisms of interactions with the other components of the environment. The following sections report a generalisation of some partial results discussed in a previous paper (Monte, 2002). 2.2. The equation controlling the population dynamics The distribution in space of biota is described by the spatial density function n(x, t, a). n(x, t, a) is the number of individuals per unit volume at point x and instant t belonging to the age-class “a” (NumberVolume−1 Age−1 ). In other words, n(x, t, a)dxda is the number of individuals that, at time t, are in a neighbourhood dx of x (here dx denotes an infinitesimal volume) and belongs to the age interval (a, a + da). Each individual traces out a path in space as time varies. Let () denote the position at instant  of a living organism that, at time t, reaches point x ((t) = x). Let us introduce the conditional probability p(x, t, a; , )d that an individual was in a neighbourhood d of  at time  assuming that it belongs to age-class “a” and is in x at time t. The biota dispersal is described by a probability function to account for the possible randomness of the movement of the organisms. In principle, the equations of the biota motion should account for the movements caused by external stimuli related to the distributions in the environment of the resources for survival and reproduction of species. We note that the use of conditional probability does not prevent us to consider the deterministic case (corresponding to probability 1). For instance, the translation with constant velocity can be represented by a probability distribution function p(x, t, a; , ) = ı(x −  − v(t − )) where ı is the Dirac’s function. p(x, t, a; , )n(x, t, a), the spatial density of individuals in x and of age-class “a” at time t that were in  at instant  (Appendix A), satisfies the following equation: ∂p(x, t, a; , )n(x, t, a) ∂p(x, t, a; , )n(x, t, a) + = (1) ∂t ∂a x (x, t, a; , ) − Mp (x, t, a; , ) + ı(a)p(x, t, a; , )LNp (x, t) where x is the movement term (net balance of organisms per unit volume around the position x, per unit age, per unit volume around the position of origin  and per unit time; Number Volume−1 Age−1 Volume−1 Time−1 ) and Mp is the biota death rate (number of organisms, per unit volume around the position x, per unit age and per unit volume around the position of origin , that die per unit time; Number Volume−1 Age−1 Volume−1 Time−1 ).

3281

LNp (x, t) is the rate of birth of the population (Number Volume−1 Time−1 ). The symbols (x, t, a; , ) indicate that the corresponding quantity refers to the individuals in x and belonging to age-class “a” (the age-class interval (a, a + da)) at time t that were in  at time . In Eq. (1) the increment of population due to the newborn individuals affects only age-class 0. Eq. (1) can be obtained by evaluating the rates of change subsequent to the movement and the population dynamics (death, reproduction and ageing) of the number of the individuals of certain age-class that were in  at time  and occupy a volume V at time t. If it is possible (Hypothesis H) to find two operators Lx and LMp (these operators are assumed to act on the space of functions of x and have codomains in the same space) independent of  and  such that, fora ≥ t −  ≥ 0 (from now on, for the sake of simplicity, the symbol d denotes the definite integral over the whole domain of the variable ):



x (x, t, a; , )d = L x n(x, t, a)

(2)

Mp (x, t, a; , )d = L Mp n(x, t, a)

(3)



by integrating Eq. (1) over the whole domain of , we get: ∂n(x, t, a) ∂n(x, t, a) + = L x n(x, t, a) − LMp n(x, t, a) + ı(a)LNp (x, t) ∂t ∂a (4) that is a generalisation of the most common form of the equation controlling the dynamics and the movement of a population (Gurtin and MacCamy, 1977, 1981; Hernández, 1998). We note that Lx and LMp may depend on x, t and a. To obtain Eq. (4) we have accounted for the normalisation condition



p(x, t, a; , )d = 1

(5)

Now we turn to the question of the existence of operators Lx and LMp . Do such operators satisfying Eqs. (2) and (3) exist? It would be helpful to find the general conditions for their existence. On the other hand, it would be very disappointing if they did not exist at all. For the purpose of the present work, it may be useful to show that they surely exist if we assume the following hypotheses that underpin several state-of-the-art ecological models: Hypothesis H1.

x and Mp can be written as follows

x (x, t, a; , ) = L x1 [p(x, t, a; , )Lx2 n(x, t, a)] Mp (x, t, a; , ) = L Mp1 [p(x, t, a; , )LMp2 n(x, t, a)]

(6)

where Lx1 , Lx2 , LMp1 and LMp2 are independent of  and  and, moreover, Lx1 and LMp1 are linear (Lx2 and LMp2 are not necessarily linear). It is trivial to demonstrate that, if Hypothesis H1 is satisfied, operators Lx and LMp exist and Eq. (4) can be derived from (1). Moreover, Lx = Lx1 Lx2 and LMp = LMp1 LMp2 in view of the normalisation condition (5) and of the linearity of Lx1 and LMp1 (Appendix B). For instance, we can write



Lx1 [p(x, t, a; , )Lx2 n(x, t, a)]d

 = Lx1

p(x, t, a; , )L x2 n(x, t, a)d = L x1 Lx2 n(x, t, a)

(7)

A similar result is valid for LMp . The biomass per unit volume is f(x, t, a) = n(x, t, a)(x, t, a) where (x, t, a) is the ensemble average (average performed over the set of individuals of age “a” in x at instant t) of the mass of the

3282

L. Monte / Ecological Modelling 220 (2009) 3280–3290

organisms. Therefore, Eq. (4) can be written as follows:





∂f (x, t, a)/ (x, t, a) ∂t = Lx 

f (x, t, a)



+

∂f (x, t, a)/ (x, t, a)

 − LMp 

(x, t, a)



∂a f (x, t, a)

 + ı(a)LNp (x, t)

(x, t, a)

(8)

Note that (x, t, a) is function of age “a” and of time t as, in principle, the mass of the organisms of a given age can be different at different instants of time due to possible environmental changes influencing the development of the individuals. In other words, the masses of individuals of a given age “a” belonging to distinct generations can be different. One of the most common examples of Lx is the dispersiontransport operator that is used when the biota movement can be approximated by a “random walk” with drift. In this case Lx2 is the identity operator and Lx1 = div(−v + Dgrad) (v and D are the drift velocity and the dispersion coefficient, respectively). Examples of LMp1 and LMp2 are multiplicative factors and “functions” of n(x, t, a), respectively. LMp and LNp can be strongly nonlinear in the population density. Indeed, the population growth, the reproduction and the mortality rates are not necessarily linearly dependent on the population size. As the number of individuals increases, the population runs out of some environmental resources. The effect on the reproduction and the mortality rates of this environmental resistance can increase with the population size (Odum, 1959). It should be noted that the population equation is frequently integrated over the domain of age “a”. Furthermore, the death, the birth and the growth terms in the equation are often aggregated and expressed as total population growth and decay and are related to the biomass density of other interacting species or environmental components (for instance, resources, preys, predators) (Dubois, 1978). Hypothesis H is essential to obtain the population equation in its most common form (Eq. (4)). It is useful to examine, briefly, the domain of application of Eq. (4) in view of some important features of the processes here considered. The terms x and Mp may depend on the “story” of the organisms or, specifically, on the paths the individuals followed to reach point x. For instance, in case of the so-called persistent random walk, an angular correlation between successive steps can be assumed (Wu et al., 2000). Furthermore, some specific features of the animal displacement, such as the leptokurtic distribution, require the application of more sophisticated mathematical approaches in the equation controlling the population dispersal (Zhang et al., 2007). It may be difficult to derive the forms of operators Lx and LMp that satisfy Eqs. (2) and (3). Moreover, in general, it is not ever possible to demonstrate that such operators exist as the first members of the above-mentioned equations might depend on the proportions of population sub-classes characterised by different “stories” (for instance, on the position that the sub-classes occupied at instant  1 ,  2 ,  3 ,  4 , etc.) rather than on the total population density n(x, t, a) in x at instant t. In view of operators Lx1 and LMp1 are linear and independent of , Eqs. (6) (Hypothesis H1) assume that the movement and death rate of the whole population, in x at instant t, are independent of the proportions of the number of individuals that at previous instants  occupied different positions  (the movement and the rates are independent of the “story” of the population, namely, of the spatial distribution of the population at instant ). In conclusion, when the movement of animals, the death and birth rates of the individuals in x at instant t depend on previously occupied positions, it is not possible, in principle, to derive Eq. (4) (Eq. (8)) in its given form.

Probabilistic models based on a large number of simulations reproducing the behaviour of the individual animals can be used in such circumstances. Before concluding this section we recall some mathematical properties of functions f and p. The above-described approach can be generalised by introducing the conditional probability p(x, t, a; 1 , 1 ; 2 , 2 ; 3 , 3 ; 4 , 4 ; . . .)

(9)

that an individual, being in x at instant t, was in  1 at instant  1 , in  2 at instant  2 ,  3 at instant  3 ,  4 at instant  4 , etc. In other words, p(x, t, a;  1 ,  1 ;  2 ,  2 ;  3 ,  3 ;  4 ,  4 ; . . .) is the probability that an organism in x at t has followed the path ( 1 ,  1 ;  2 ,  2 ;  3 ,  3 ;  4 ,  4 ; . . .). An equation similar to (1) can be written by using the conditional probability for each given path. Eq. (4) can be derived provided that the previous operations of integration are performed over all the variables  1 ,  2 ,  3 ,  4 , etc. n(x, t, a) tends to 0 as “a” tends to infinity (or, in terms of a more strict condition, n(x, t, a) = 0 if a > M where M is a sufficiently large value). This assumption states that the lifetime of any individual is finite. The use of the Dirac’s delta function in the birth rate term requires the application of the notion of generalised functions of variable “a”. n(x, t, a) shows a discontinuity at a = 0 (the function jumps from 0 to a value n0 ) and satisfies the condition n(x, t, a) = 0 when a < 0. Consequently,



∂ n(x, t, a)da = 0 ∂a Da

(10)

where Da is the whole domain of the age-class. For any function F of  and , we can write:



p(x, t, a; 1 , 1 ; 2 , 2 ; 3 , 3 ...; , )F(, )dd1 d2 d3 . . .

 =

(11) p(x, t, a; , )F(, )d

The above mathematical results will be used in some applications that will be illustrated in the following sections. 2.3. Modelling the behaviour of contaminant concentrations in migrating biota We will derive the equation controlling the contaminant migration by following each individual along its trajectory and by averaging the contaminant concentration over the ensemble of individuals that have reached the generic point x (Lagrangian method, Appendix C). Let us introduce the vector q(x, t, a; ) whose components are the amount of contaminant in different organs at time  in a particular individual that moves along a given trajectory () and that at time t is in x and belongs to age-class “a” (more specifically, each component qi of the vector is the mass of contaminant in a particular organ or tissue–liver, bone, blood, etc.). We assume that q(x, t, a; ) obeys the following equation: dq(x, t, a; ) = B q(x, t, a; ) + ℘(, , a − t + ) d

(12)

℘(, , a − t + ) is a vector whose generic component j is the mass of pollutant migrating, per unit time, at instant  from the external environment (including abiotic and biotic components) to tissue (organ) j of the individual in x and belonging to age-class “a” at instant t (therefore, the individual belongs to age-class a − t +  at instant ). B is an operator acting on the space of vectors q (the jth component of B q(x, t, a; ) is the net rate of variation of contaminant in organ j due to exchanges with other organs and to the

L. Monte / Ecological Modelling 220 (2009) 3280–3290

3283

Fig. 1. Pictorial representation of trajectories. The subscripts b, c, d . . . denote individuals.  b () is the position at time  of an individual of age “a” in x at instant t. qb (x, t, a; ) is the amount of contaminant at instant  in the individual “b” that at time t is in x and belongs to age-class “a”. b (a − t + , ) is the mass of the individual at instant . The ratio (qb (x, t, a; ) + qc (x, t, a; ) + qd (x, t, a; ) + qe (x, t, a; ) + ···)/(b (a − t + , ) + c (a − t + , ) + d (a − t + , ) + e (a − t + , ) + ···) = /<> at instant t ( = t) is C(x, t, a) (average concentration).

release to the external environment or due to any kind of chemical/physical degradation process such as the radioactive decay). As usual,  is the position of the individual at time . B is a stochastic operator and ℘(, , a − t + ) is a stochastic vector as they depend on the specific individual. It is important to note some features of Eq. (12). ℘(, , a − t + ) is commonly related, by means of suitable sub-models, to the average concentrations of pollutant in the environmental components interacting with the individual: ℘(, , a − t + ) =



(,,a−t+,˛)

Ki

Wi (, , ˛)d˛

(13)

i

where Ki are operators that act on functions of . Wi (, , ˛) is the spatial concentration of pollutant in the ith environmental component in  at instant  belonging to age-class ˛. Functions Wi (, , ˛) depend on the age-class ˛ as the external components can be other living organisms. The superscript symbols (,,a−t+,˛) emphasise that the processes of incorpoin operators K i ration of pollutants in an organism depend on the spatial position (), the time (), the age-class (a − t + ) of the organism at instant  and on the age-class ˛ of the other environmental components (for instance the age of preys in case of selective predation depending on the prey size). A trivial example of application of Eq. (12) is the compartment model. In that case, the operator B is the matrix of pollutant migration constant rates between organs and from the whole organism to the external environment and Ki are multiplicative factors related to the coefficients of migration from the external environmental components to the tissues of the individual organism (Rescigno et al., 1983). Fig. 1 depicts the trajectories and the quantities previously defined. To manage Eq. (12) we hypothesise that it is possible to find an (,,a−t+) operator L B and a vector M(, , a − t + ) such that: Hypothesis H2.







B q(x, t, a; ) =

and



℘(, , a − t + ) =

(,,a−t+)

p(x, t, a; , )L B

Q (, , a − t + )d (14)

 p(x, t, a; , )M(, , a − t + )d

(15)

where Q(x, t, a) is the ensemble average of q(x, t, a; t) over the set of animals that at time t are in an “infinitesimal” neighbourhood of

x and belong to age-class a: Q(x, t, a) = q(x, t, a, t)

(16)

In other words, Q(x, t, a) is the average amount of contaminant in the animals that at time t have reached the point x whatever the trajectory () was. From Hypotheses H1 and H2, it is possible to demonstrate, by boring calculations (Appendix C), that the contaminant obeys the following equation: ∂Q (x, t, a)n(x, t, a) ∂Q (x, t, a)n(x, t, a) + = n(x, t, a)L B Q (x, t, a)+ ∂t ∂a +Lx1 [Q (x, t, a)Lx2 n(x, t, a)] − LMp1 [Q (x, t, a)LMp2 n(x, t, a)]+ +ı(a)Q (x, t, a)LNp (x, t) + n(x, t, a)M(x, t, a) (17) Eq. (17) is a generalisation of the model to predict the dynamics of the population–toxicant interaction (Buonomo et al., 1999). To put Q(x, t, a) in place of p(x, t, a; , ) in Eqs. (1) and (6) is the rule of thumb to include the population dynamics terms in Eq. (17). As the concentration of contaminant can be defined as (x,t,a) C(x, t, a) = Q , Eq. (17) can be written in terms of C(x, t, a) (x,t,a)

and f(x, t, a). Vector M and operator LB surely exist if, as commonly hypothesised by several models (Appendix D), we can assume that:

Hypothesis SH1. ℘(, , a − t + ) depends on position  and ageclass at instant  but is independent of the individual animal (it is not stochastic); therefore M(, , a − t + ) = ℘(, , a − t + ) Hypothesis SH2a. B is linear and can be assumed independent of the individual animal, or, alternatively Hypothesis SH2b. B is independent of the individual animal (it should be not necessarily linear) and for each animal in x at instant t q(x, t, a; ) = Q (, , a − t + ), where  is the position of the animal at time  (the amount of contaminant, at instant , in an individual that is in x at time t and has followed the trajectory () is exactly equal to the “ensemble” average at position  an time ) (Appendix D). If Hypothesis SH2a or Hypothesis SH2b are verified LB = B . It is important to note that Hypotheses SH2a requires, moreover, that the ensemble average of contaminant content q at time  over sub-classes of individuals of age “a” in x at instant t that were in  at time  is equal to the average over the entire class of individuals in  at instant  and of age (a − t + ) (therefore, regardless the positions at instant t).

3284

L. Monte / Ecological Modelling 220 (2009) 3280–3290

Although the mathematical details are tedious, may seem difficult to understand and too demanding to read, the meaning of the above-mentioned hypotheses can be easily summarised and interpreted in an ecological perspective. Hypothesis H and, consequently, Hypothesis H1 state that the movement and the death rate of individuals at a given instant t can depend on the position x and the age-class “a” at that instant, but are independent of the positions previously reached by the organisms (more generally they are independent of the individual “story”). Hypothesis SH1 states that the rates of uptake of contaminant from the environmental components, in a given point and at a given instant, must be approximately equal for all the individuals in the same age-class. Similarly, Hypotheses SH2a and SH2b state that the dynamics of contaminant in different individuals belonging to a given age-class must be, approximately, the same (the rates of migration between organs and from the individual to the external environment should be approximately equal for different individuals). In other words, Hypotheses SH1, SH2a and SH2b assume that the effects of the biological variability among individuals on the uptake and excretion mechanisms are approximately negligible. Hypothesis SH2b states, moreover, that the concentration in an individual that at a given instant occupies a given position and belongs to a given age-class is approximately equal to the average concentration calculated over all the animals around the same position at same time and in the same age-class. In other words, it should be possible to assume that the variability of contaminant concentration among individuals in a given age-class and around a certain position at a given instant is, approximately, negligible. Note that the linearity in Hypothesis SH2a is valid when the biological processes of contaminant migration are simulated by the commonly used compartment models. In conclusions, Hypotheses H1 and H2 allow one to write Eq. (17). Hypothesis H2 is certainly valid if Hypotheses SH1 and SH2a or SH2b are true. Hypotheses H1, SH1 and SH2a or SH2b assert that sub-classes of operators Lx , LMp , LB and of vectors M, for which Hypotheses H1 and H2 are satisfied, surely exist. Luckily, Hypotheses H1, SH1 and SH2a or SH2b underpin several state-of-the-art models for predicting the behaviour of population dynamics and of contaminant in abiotic and biotic environmental components. It is easy to realise that the Hypotheses H1, SH1 and SH2b are assumed valid for continuous or quasi-continuous media like water and particulate matter suspended in water (the heuristic assumption is that the “individuals” transporting contaminant are infinitesimal volumes of media). 2.4. Coupling the equations If we combine Eqs. (4) and (17) with the sub-models for assessing the birth and the death rates as functions of the environmental contamination and/or of the internal contamination of organisms, we obtain a comprehensive model for predicting the dynamic of biota populations following the detrimental effects of environmental pollution. It is quite obvious that the procedure can be very complicated and that a large amount of data and information, which are not ever easy to obtain, are necessary. Moreover, as we have noted in the previous sections, it is not possible to predict the contaminant behaviour in biota and the biota population dynamics by Eqs. (4) and (17) when the movement of organisms at a given instant is related to the past positions. As we have previously remarked, in such circumstances it is necessary, at least in principle, to make use of different kinds of predictive tools such as stochastic models that, by means of a sufficient number of runs, simulate the movement and the contamination dynamics of a great deal of different animals.

In case of radioactive contamination a simplified strategy based on effect–dose assessments was suggested. The approach assumes causal relationships between the values of quantities measuring the risk associated with ionising radiation, such as the dose rates, and the parameters in the equations controlling the dynamics of the populations of aggregated ecosystem components. Examples of application to environmental areas highly contaminated by radionuclides (Monte, 2009a) showed that this kind of models can explain certain systemic responses of ecosystems polluted with radioactive substances.

3. Discussion 3.1. Domain of application of the methodology Generally, the migration of a contaminant through the environment is the result of the transport by more or less mobile carriers. Obvious examples of carriers are masses of air or water, particulate matter and biota. Firm substrates, like bottom sediments or soils, are further examples of “carriers” whose translation velocity is usually assumed to be zero. A carrier is composed of elemental units, the “carrier elements” which move through the environment. The methodology presented in the previous sections can be applied also in the case the “carriers” transporting the pollutant are not necessarily living organisms. It is sufficient to change the words “animal” or “living organism” with the term “carrier element” in any sentence. The “age-class” is commonly used to characterise the age distribution in a population of living organisms. Nevertheless, it can be conveniently applied to non-living objects. The concept of “ageclass” can be extended to abiotic components assuming that the age is the length of the existence of a carrier element extending from the moment the element entered the environmental system to instant t. For instance, we can attribute an age-class to the inlet water inflowing into a lake by assuming that the age “a” of a water element at instant t is the time extending from the moment the element entered the lake to time t. We can assume that the overall dispersion of a contaminant in the environment is due to the movement of the different kinds of carriers that transport the contaminant itself. As previously stated, the terms birth rate and death rate have a broader meaning than the corresponding biological notions. Indeed, they should be intended as the rates at which masses of carrier elements enter and leave the modelled system per unit volume (Eq. (8)). The amount of contaminant in an environmental component is generally measured as the average over the carrier elements sampled in a finite region of space. In other words, the measure of the amount of a pollutant is a “macroscopic” quantity corresponding to the ensemble average Q(x, t, a) (Eq. (16)). The amount of contaminant in a single carrier element is the “microscopic” value. For instance, the “macroscopic” quantity of a contaminant in particulate matter suspended in water or in atmosphere is the total amount of contaminant in a sample of particulate matter divided by the number of the particles in the sample. The “microscopic” quantity, q(x, t, a; ), is the amount of pollutant in each single particle. Similar definitions can be applied to the macroscopic and the microscopic concentrations. To speak of macroscopic and microscopic quantities or concentrations depicts realistically the case for particulate matter and for many other types of carriers of small size. Difficulties may arise due to the impossibility of defining appropriately the ensemble average in some situations. This can be indeed meaningless in case of few large moving living organisms.

L. Monte / Ecological Modelling 220 (2009) 3280–3290

In such circumstances it is necessary to follow each single animal along its trajectory. It is worthwhile to note that Hypothesis SH2b states that the “macroscopic” amount of pollutant at a given point is equal to the corresponding “microscopic” quantity at any time t (q(x, t, a; t) = Q(x, t, a) where we put  = t). 3.2. Comparison with mass balance approach The derivation of Eq. (17) by the mass balance method is commonly used in case of abiotic carriers. It is normally based on the assumption that the dynamics of the contaminant in the carrier elements, the carrier movement, the rates of contaminant exchange between the carrier elements and the external environmental at a given position x and instant t depends on x and t but are independent of the “story” of the contaminated carrier element. The above hypotheses are commonly accepted for abiotic matter. Indeed, Eq. (17) is the prototype of the most common models used to predict the behaviour of contaminant in non-living components of the environment (Monte, 2009b). The analysis performed by applying the Lagrangian method better enlightened that Eqs. (4) and (17) cannot be generally applied in the case of living organisms as the population dynamics, the biota movements and the dynamics of pollutant can be significantly influenced by the “story” of the individuals. 3.2.1. “Ambiguities” in the movement term in Eq. (17) To assume that the dispersion of a carrier is controlled by a Fickian diffusion mechanisms is common to several models. It is not the aim of the present work to discuss in which circumstances the approach is valid. We will rather try to discuss some ambiguities that the mass balance approach, at first sight, does not clearly enlighten and that the present methodology helps to explain and remove. It is apparent that, assuming that both Lx1 and Lx2 are linear, from the following two different forms of the movement term x (x, t, a; , ) = L x1 [p(x, t, a; , )Lx2 n(x, t, a)]

(18a)

x (x, t, a; , ) = L x1 Lx2 [p(x, t, a; , )n(x, t, a)]

(18b)

by an integration over the  domain, we derive, for both expressions (18a) and (18b), the same equation controlling the population dynamics (Eq. (5)). On the contrary, the correspondent terms that control the dispersion of the contaminant transported by the carrier are different in Eq. (17) as they depend on the particular composition of operators Lx1 and Lx2 in (18a) and (18b): L x1 [Q (x, t, a)Lx2 n(x, t, a)]

(19a)

L x1 L x2 [Q (x, t, a)n(x, t, a)]

(19b)

In particular, in the case of Fickian diffusion mechanisms, two different forms of the term for calculating the pollutant flux due to the dispersion are found in the scientific literature (D denotes the dispersion coefficient): Buonomo et al. (1999): div[Q (x, t, a)Dgradn(x, t, a)]

(20a)

Zheleznyak et al. (1992): divDgrad[Q (x, t, a)n(x, t, a)]

(20b)

If (x, t, a) is independent of the position x, in view of the (x,t,a) linearity of operator grad and accounting for C(x, t, a) = Q (x,t,a)





and f (x, t, a) = n(x, t, a) (x, t, a) , it is possible to write div[C(x, t, a)Dgradf (x, t, a)]

(20c)

divDgrad[C(x, t, a)f (x, t, a)]

3285

(20d)

Eqs. (20c) and (20d) model the migration of a pollutant for two different kinds of motion of the carriers. The first equation is valid when the movement of the sub-class of individuals that were in a given position at a given instant depends on the density of the whole population in x, whereas, the second equation can be applied when the movement depends on the density of the population of the subclass. To enlighten the consequences of this important difference, we can consider the simple case of a population density with null gradient. The flux of contaminant −C(x, t, a)Dgrad f(x, t, a) is 0 in the first case (20c), whereas it can be different from 0 in the second case (−Dgrad[C(x, t, a)f(x, t, a)]) if the concentration of pollutant is not constant with x. For non-living carriers subject to Fickian dispersion, it is commonly assumed that expression (20d) is the appropriate one (Zheleznyak et al., 1992). On the contrary, for living organisms, it is possible that the movement of an individual is controlled by the whole population density. Consequently, if the population is uniformly distributed (n(x, t, a) = constant), from Eq. (19a), assuming Lx2 = grad, it follows that the movement term is null in Eq. (1) which controls the behaviour of p(x, t, a; , ) n(x, t, a). Therefore, the displacement of each single individual is confined within a neighbourhood of a given position x. Eq. (18a) better models the movement of species with a marked territorial behaviour. On the contrary, (18b) assumes that each single individual can move freely throughout the area occupied by the whole population. 3.3. Examples of applications We have seen that the use of Eqs. (4) and (17) shows several limitations in view of the many hypotheses we have assumed. It is therefore useful to illustrate examples of possible successful applications of these equations. 3.3.1. Contamination of moving biota Commonly, state-of-the-art models that include the effects of animal movement (worms, insects, etc.) on the contaminant dispersion in soil or sediment (“bioturbation”) assume that this process produces an enhanced dispersion of contaminant. Eq. (17) offers the opportunity for approaching this problem within a more general perspective. We will describe, briefly, an example of application of the method to the contamination of salmons undertaking spawning migrations. This kind of problem can be easily approached by solving the equations controlling the contaminant behaviour in a frame of reference moving with the fishes. Therefore, this particular application can be a suitable example to check the methodology here described. We will assume that only the class of adult salmons is involved in the population migration. Therefore we will disregard the dependence on the age-class of functions in Eq. (17). For simplicity we will consider the one-dimensional case (x is the distance of the salmon school from the river mouth). During the migration the population density of salmons in the watercourse can be represented by a function f(x, t) that is the biomass of salmons in x at time t per unit length. If R is the mortality rate (Time−1 ) during the migration, the population dynamics (Eq. (8)) is described by the following equation (we have assumed that the mass of salmons is, approximately independent of x and t): ∂f (x, t) ∂f (x, t) = −v − Rf (x, t) ∂t ∂x

(21)

3286

L. Monte / Ecological Modelling 220 (2009) 3280–3290

where we have used the translation operator Lx = −v

∂ ∂x

(22)

to model the fish movement. The general solution of Eq. (21) is f (x, t) = f ∗ (x − vt)e−Rt

(23)

where f* is an arbitrary function representing the initial population density. To illustrate this simple application, we assume that the complex processes of contaminant exchange between the water and the fish can be approximately modelled by a first order differential equation: ∂q(x, t) = − B q(x, t) + K ∗ W(x, t) ∂t

(24)

where q(x, t) is the amount of contaminant in the single individual, W(x, t) is the contaminant concentration in water, B is the contaminant excretion rate constant (s−1 ) and K* is the contaminant transfer factor from water to biota (m3 s−1 ). This is the simplest and the most frequently used model for predicting the behaviour of a contaminant in fish species. By writing Eq. (17) in terms of f(x, t) and C(x, t) (the concentration of contaminant in salmons), we obtain: ∂C(x, t)f (x, t) ∂C(x, t)f (x, t) = − B f (x, t)C(x, t) − v + ∂t ∂x −RC(x, t)f (x, t) + f (x, t)KW (x, t)

(25)

where K = K* / ( is the average mass of a fish). Eq. (25) states that the spatial concentration (C(x, t)f(x, t), mass per unit volume) of contaminant transported by the salmon school in the river is subject to a translation movement with velocity v. The methodology can be also applied to other processes of environmental contamination that involve the movement of animals from more or less contaminated areas. This can occurs, for instance, in deep lakes when surface waters (epilimnion) and deepwater layers (hypolimnion) show different levels of contamination (Monte et al., 1993). 3.3.2. Age-class effect The metabolism of a contaminant depends on the age of the living organism. Alimentary habits of animals can change during the life depending on the stage of development of the individuals. Moreover, the concentration of contaminant in biota depends on the different periods of exposure to the contamination of individuals born at different times. An increase of the variability of the contamination among individuals following an accidental contamination of the environment and a linear correlation between the contamination and the ageclass of fishes at a given time were observed in perch in Lake IJsselmeer (The Netherlands) contaminated by 137 Cs introduced into the environment following the Chernobyl accident (Heling, 1997; IAEA, 2000). To focus the analysis specifically on the effect that different times of exposures to the contaminant may have on the contamination levels in biota, we will use, for simplicity, the biological model (24). The population equation is as follows: ∂f (t, a) ∂f (t, a) + = −Rf (t, a) ∂t ∂a

(26)

where R is the death rate of the fishes. Let us assume that (a) the radionuclide is homogeneously distributed in the lake (consequently we will not account for the movement term in Eq. (17));

Table 1 Slope and intercept of the linear relation between the concentration of 137 Cs and the age of perches sampled in Lake IJsselmeer following the fall-out from the Chernobyl accident (IAEA, 2000). The first column reports the year of sampling. Year

Km (Bq kg−1 year−1 dry weight)

Kp (Bq kg−1 dry weight)

Correlation coefficient r2 between contamination and year-class)

1988 1989 1990

13.1 9.1 5.1

12.7 6.7 6.4

0.837 0.586 0.661

(b) the radionuclide concentration in water, CW , can be approximately fitted to an exponential decay curve CW (t) = W0 e−At in the period 1988–1990 (IAEA, 2000); (c) the initial amount of contaminant in fish fry at time t is approximately proportional to the concentration in water, consequently C(t, 0) = C0 e−At . Eq. (17) becomes: ∂C(t, a)f (t, a) ∂C(t, a)f (t, a) + ∂t ∂a = − B C(t, a)f (t, a) − RC(t, a)f (t, a) + f (t, a)KCW (t)

(27)

The solution of Eq. (27), when a < t, is: C(t, a) =

KW0 e−At (1 − e−( B −A)a ) + C0 e−A(t−a) e− B a B − A

(28)

If the product ( B − A)·a is small, Eq. (28) becomes: C(t, a) = KW0 e−At a + C0 e−At

(29)

The half-life of the concentration of contaminant in water and in perch were both of the order of 1 year in Lake IJsselmeer (IAEA, 2000). Therefore it is reasonable to hypothesise that formula (29) is approximately valid. The analysis of the available empirical data demonstrated that the concentrations Cp of 137 Cs in perch in Lake IJsselmeer sampled in years 1988, 1989 and 1990 can be approximately fitted to a linear function of the age of the fish: Cp = Km a + Kp where Km and Kp are two empirical parameters (IAEA, 2000). Table 1 shows the values of the parameters and of the correlation coefficient (r2 ) calculated for fishes sampled in different years. Eq. (29) states that the concentration is a linear function of the age and that the angular coefficient KW0 e−At and the intercept C0 e−At of the corresponding line decrease with time, in agreement with the empirical data. 3.3.3. A counterexample We will show now a counterexample that illustrates one of the possible shortcomings of Eq. (17). Eq. (17) contains a first order derivative with respect to time. This reflects the expected result that the knowledge of the initial spatial distribution of the contaminant is a sufficient condition for predicting the time behaviour of the contamination when the particular mechanisms of migration are specified (that is, the form of the terms in the right side of the equation). As previously stated, the “macroscopic” (ensemble average) and the “microscopic” quantities of contaminant in a carrier are not necessarily equal. For instance, we noted that the experimental measure of contaminant concentration in an environmental component is a “macroscopic” quantity for it is obtained by sampling the carrier over a finite region of space. It is obvious that different distributions of contaminant at microscopic level can correspond to the same particular value of the measured “macroscopic” quantity. For instance, a certain value of the macroscopic contamination, say an average of 100 mg individual−1 , is measured when:

L. Monte / Ecological Modelling 220 (2009) 3280–3290

(a) 50% of carrier elements show a “microscopic” contamination of 200 mg individual−1 and the contaminant in the remaining 50% is negligible; (b) the contaminant in 100% of the elements is 100 mg individual−1 . Let us assume that: 1. the carrier elements are sufficiently randomly distributed so that collected sample shows an average concentration of 100 mg individual−1 ; 2. B is non-linear (for instance B q = − q2 ). In that case, neither Hypothesis SH2a nor Hypothesis SH2b in Section 2.3 are verified. For simplicity, we assume, moreover, that B is independent of the specific carrier element (in other words it is not stochastic). The variation of the macroscopic contamination at instant t is −2002 dt (n/2)/n, in case (a), and −1002 ndt/n, in case (b) (n is the number of sampled individuals). In conclusion, although the macroscopic average amount per individual at instant t are equal, their values at instant t + dt are different. Consequently, the pollutant migration process cannot be predicted by Eq. (17) that, being of the first order, has a unique solution when the initial spatial distribution of pollutant is assigned. 4. Conclusions The scope of this work was the presentation and the discussion of a methodology rather than of specific environmental models. The described general conceptual approach for predicting the migration of contaminant substances through the environment assumes that the contaminant is transported by certain mobile carriers (both biotic or abiotic). By applying the particle-tracking methodology, it is possible to derive Eqs. (4) and (17) that are the basis of several state-of-theart models for predicting the behaviour of pollutants transported through the environment by abiotic (such as particulate matter, water, air) (Monte, 2009b) or biotic carriers. We have demonstrated, by some examples, that the described methodology, in spite of the limitations of its domain of application, can be used in several circumstances. Moreover, the methodology offered the opportunity of settling an ambiguity in the dispersion term of the pollutant migration equation (Section 3.2.1). The derivation of the equations by the Lagrangian method enlightens the hypotheses necessary for the proper application of the methodology here presented. These hypotheses are not readily apparent when the equations are derived by the mass balance approach. The principal assumption is that the movement of an individual, at a given instant, must not depend on the “story” of the individual itself. This is commonly assumed for abiotic matter although it is not ever true for biota. Consequently, it is not ever possible to model, by Eqs. (4) and (17), the spatial dispersion of a population of living organisms and the dynamics of the contaminant concentration in these organisms when the dynamics of movement at a certain instant is a function of the previous conditions. Furthermore, we have shown that the characteristics of the mechanisms of interaction of the pollutant with the carrier elements (chiefly the non-linearity with pollutant concentrations and the variability among individuals (Section 3.3.3)) do not ever allow to derive Eq. (17) even when Eq. (4) can be assumed valid. A common alternative approach to overcome the above difficulties is based on the application of probabilistic techniques, such as

3287

Monte-Carlo methods, to derive the spatial distributions of populations and pollutants by simulating the behaviour of a large number of individuals. Unfortunately, such a method can be very computing time-expensive. In general, the application of models developed by the above mentioned methodologies requires the knowledge of a great deal of data and information. As frequently noted, simplified approaches based on process aggregation can be a further strategy to approach this kind of complex environmental modelling problems. Appendix A. Let n(x, t, a)dx denote the number of individuals of age “a” in a neighbourhood of x at instant t. n(x, t, a; , )dxd, the number of individuals in x at time t that were in  at instant , is n(x, t, a; , )dxd = p(x, t, a; , )n(x, t, a)dxd

(A1)

The biomass spatial density f(x, t, a) of individuals in x at instant t is:





f (x, t, a) = n(x, t, a) (x, t, a)

(A2)

where (x, t, a) is the average mass of individuals in x at instant t. It is important to note that all the equations in this work can be f (x,t,a) written in terms of f(x, t, a) by substituting n(x, t, a) with (x,t,a) . Appendix B. We briefly recall the notion of linear operator. An operator acting on functions f(x) is linear if L(f (x) + g(x)) = Lf (x) + Lg(x) L(cf (x)) = cL(f (x))

(1B)

where c is a number and f(x) and g(x) are two functions. If L is linear and independent of  and  we can write

 



Lf (x, , )d = L Lf (x, , )d = L



f (x, , )d (2B) f (x, , )d

(2B) is valid under some hypotheses that are generally true for the kinds of functions and operators commonly used in environmental modelling. The transport-diffusion operator is a typical example of linear operator. A similar definition of linearity can be applied to operators acting on vector spaces (B ). Appendix C. (,,a−t+)

In Section 2.3, we wrote L B to recall that, in principle, LB may depend on time  on the age-class (a − t + ) and on the position  of the carrier element at instant . The occurrence that operator LB and function ℘(, , a − t + ) may depend on the age classes of the carrier element is not surprising. For instance, the metabolism of biota is a function of the individual age. LB and ℘(, , a − t + ) can explicitly depend on  as the complex interactions of the pollutant with the carrier may be influenced by several specific environmental parameters and conditions (temperature, chemical factors, etc.) that can vary with time. By averaging q(x, t, a; ) over all the trajectories that, at instant t, reach x (Fig. 1), Eq. (12) becomes:



∂ q(x, t, a; ) ∂





 



= B q(x, t, a; ) + ℘(, , a − t + )

(1C)

3288

L. Monte / Ecological Modelling 220 (2009) 3280–3290

We can derive the next formulae also for the more general case that probability refers to a multi-step path (p(x, t, a; ␰1 ,  1 ;  2 ,  2 ;  3 ,  3 ;  4 ,  4 ; . . .)) in view of Eq. (11). Formally, the general integral of Eq. (1C), in view of Eqs. (14) and (15) (Hypothesis H2), is:





q(x, t, a;   )



Note that, if LB depends on  and on the age-class (a − t + ), it depends on t and a in term C (Eq. (5C)) and, therefore, in Eq. (10C). Similarly, in virtue of the occurrence of p(x, t, a; .t) in the integrand of term C, LB in (10C) depends on x. In Eq. (10C), ı(a)Q (x, t, a)LNp (x, t) should be 0 as Eq. (3C) assumes that the contamination is 0 when a = 0 (Q(x, t, 0) = 0).





=

(2C)

p(x, t, a; , )(LB Q (, , a − t + ) + M(, , a − t + ))dd t−a

that at   = t gives:



Q (x, t, a) = q(x, t, a; t)









=

(3C)

p(x, t, a; , )[(L B Q (, , a − t + ) + M(, , a − t + ))]dd t−a

Multiplying both members of Eq. (3C) by n(x, t, a) we obtain: ∂n(x, t, a)Q (x, t, a) =A+B+C −D ∂t

(4C)

where



 

t

A=



t−a



t

B=

The contamination of newborns can be included in the right side of Eq. (1C) by adding the term 2q0 (x, t, a; )ı( − t + a)

(11C)

where q0 (x, t, a; t − a) is the amount of pollutant at instant t − a in a carrier element that is in x and belongs to the age-class “a” at time



∂ n(x, t, a)p(x, t, a; , ) [LB Q (, , a − t + ) + M(, , a − t + ))]dd ∂t

n(x, t, a)p(x, t, a; , )

t−a

∂ [LB Q (, , a − t + ) + M(, , a − t + ))]dd ∂t

(5C)

n(x, t, a)p(x, t, a; , t)[LB Q (, t, a) + M(, t, a)]d = n(x, t, a)[LB Q (x, t, a) + M(x, t, a)]

C=

 D=

n(x, t, a)p(x, t, a; , t − a)[LB Q (, t − a, 0) + M(, t − a, 0)]d

We have used p(x, t, a; , t) = ı(x − ) to calculate C in Eq. (5C). Similarly, ∂n(x, t, a)Q (x, t, a) = A1 + B1 + D1 ∂a





t

A1 =



t−a t





∂ n(x, t, a)p(x, t, a; , ) [LB Q (, , a − t + ) + M(, , a − t + ))]dd ∂a

B1 =

n(x, t, a)p(x, t, a; , ) t−a

 D1 =

t and occupied position  at time t − a (t − a is the instant of birth or the instant the element entered the system).

∂ [LB Q (, , a − t + ) + M(, , a − t + ))]dd ∂a

(6C)

n(x, t, a)p(x, t, a; , t − a)[LB Q (, t − a, 0) + M(, t − a, 0)]d

In virtue of [LB Q (, , a − t + ) + M(, , a − t + ))]

(7C)

is a function of a − t + , we have B + B1 = 0. If Hypothesis H1 is valid, from Eqs. (1), (5C) and (6C), it follows: A + A1 = Lx1 [Q (x, t, a)Lx2 n(x, t, a)] − LMp1 [Q (x, t, a)LMp2 n(x, t, a)] +ı(a)Q (x, t, a)LNp (x, t)

(8C)

Moreover −D + D1 = 0

Averaging (11C) over , we get:



 =2



2q0 (x, t, a; )ı( − t + a)

p(x, t, a; , )Q 0 (, ,  − t + a)ı( − t + a)d

(12C)

Q0 (, t − a, 0) is the “macroscopic” amount at time t − a of contaminant in the carrier elements of age-class 0 in  (we have used the hypothesis assumed for Eqs. (14) and (15)). After some calcu-

(9C)

Finally we obtain: ∂Q (x, t, a)n(x, t, a) ∂Q (x, t, a)n(x, t, a) + = n(x, t, a)L B Q (x, t, a)+ ∂t ∂a +Lx1 [Q (x, t, a)Lx2 n(x, t, a)] − LMp1 [Q (x, t, a)LMp2 n(x, t, a)] + ı(a)Q (x, t, a)LNp (x, t)+ +n(x, t, a)M(x, t, a)

(10C)

L. Monte / Ecological Modelling 220 (2009) 3280–3290

as p(x, t, a; t) = ı(x − ). Eqs. (7D) show that we can interpret LB as a sort of “average” of operator B acting over the average contamination Q(x, t, a) and M as the “average” of ℘(x, t, a).

lations we obtain: Q (x, t, a)

 t

p(x, t, a; , t)(LB Q (, , a − t + ) + M(, , a − t + ))dd+

=

References

 t−a +

p(x, t, a; ,  − a)Q 0 (, t − a, 0)d

(13C) It is trivial to verify, by repeating the calculations from (4C) to (9C), that (13C) satisfies (10C) and that the term Q(x, t, a)LNp (x, t) should not be necessarily null when a = 0 (it is Q0 (x, t, 0)LNp (x, t)) Appendix D. Comments to mathematical equations Let qn,m denote the amount of contaminant at instant  in an individual identified by the couple n, m and that is in x at time t. The individual n, m occupies position  m at instant . If Hypothesis SH2a is verified we have





B q(x, t, a; ) B q1,1 + B q2,1 + B q3,1 . . . + B q1,2 + B q2,2 + l B q3,2 . . . + . . . B q1,k + B q2,k + B q3,k . . . = N

(1D)

where N is the total number of animals in x and of age “a” at time t. Since B is linear, we can write





B q(x, t, a; ) B (q1,1 + q2,1 + q3,1 . . .) + B (q1,2 + q2,2 + q3,2 . . .) + . . . B (q1,k + q2,k + q3,k . . .) = N

(2D)

By using the definitions of p(x, t, a; , t) and of Q and accounting for Hypothesis SH2a we obtain





B q(x, t, a; ) = B p(x, t, a; 1 , )Q(1 , , a − t + )d + B p(x, t, a; 2 , )Q(2 , , a −t + )d + . . . . . . + B p(x, t, a; k , )Q(k , , a − t + )d + . . . (3D) or







B q(x, t, a; ) =



p(x, t, a; , )B Q (, , a − t + )d

(4D)

If Hypothesis SH2b is verified we can write:



B q(x, t, a; ) B Q(1 , , a − t + )N1 + B Q(2 , , a − t + )N2 + . . . + B Q(k , , a − t + )Nk = N

(5D)

where N1 , N2 , . . . Nk are the numbers of animals in x at instant t that were, at instant , in  1 ,  2 , . . .  k , respectively. From (5D) we can derive again (4D). Accounting for Hypothesis SH1, we can write:



3289



℘(, , a − t + ) ℘(1 , , a − t + )N1 + ℘(2 , , a − t + )N2 + . . . + ℘(k , , a − t + )Nk = N

(6D)

From (6D) we can derive (15). The meaning of Eqs. (14) and (15) can be better understood if we evaluate the averaged values at  = t

  (x,t,a) B c(x, t, a;  t) = LB Q (x, t, a) ℘(x, t, a) = M(x, t, a)

(7D)

Beresford, N.A., Balonov, M., Beaugelin-Seiller, K., Brown, J., Copplestone, D., Hingston, J.L., Horyna, J., Hosseini, A., Howard, B.J., Kamboj, S., Nedveckaite, T., Olyslaegers, G., Sazykina, T., Vives, I., Batlle, J., Yankovich, T.L., Yu, C., 2008. An international comparison of models and approaches for the estimation of the radiological exposure of non-human biota. Applied Radiation and Isotopes 66, 1745–1749. Buonomo, B., Di Lillo, A., Sgura, I., 1999. A diffusive-convective model for the dynamics of population–toxicant interactions: some analytical and numerical results. Mathematical Biosciences 157, 37–64. Crist, T.O., Wiens, J.A., 1995. Individual movements and estimation of population size in darkling beetles (Coleoptera: Tenebrionidae). Journal of Animal Ecology 64, 733–746. Dubois, D.M., 1978. State-of-the-art of predator-prey system modelling. In: Jørgensen, S.E. (Ed.), State-of-the-art in Ecological Modelling, Proceeding of the Conference on Ecological Modelling. Copenhagen, Denmark, pp. 163–217. Goodwin, R.A., Nestler 1, J.M., Anderson, J.J., Weber, J.L., Loucks, D.P., 2006. Forecasting 3-D fish movement behavior using a Eulerian–Lagrangian-agent method (ELAM). Ecological Modelling 192, 197–223. Gurtin, M.E., MacCamy, R.C., 1977. On the diffusion of biological populations. Mathematical Biosciences 33, 35–49. Gurtin, M.E., MacCamy, R.C., 1981. Diffusion models for age-structured populations. Mathematical Biosciences 54, 49–59. Håkanson, L., Bryhn, A.C., 2008. A dynamic mass-balance model for phosphorus in lakes with a focus on criteria for applicability and boundary conditions. Water Air Soil Pollution 187, 119–147. Håkanson, L., Peters, R.H., 1995. Predictive Limnology. Methods for Predictive Modelling. SPB Academic Publishing, Amsterdam, p. 464. Hallam, T.G., Clark, C.E., Jordan, G.S., 1983. Effects of toxicants on populations: a qualitative approach II. First order kinetics. Journal of Mathematical Biology 18, 25–37. Heling, R., 1997. LAKECO: modelling the transfer of radionuclides in a lake ecosystem. Radiation Protection Dosimetry 73, 191–194. Hernández, G.E., 1998. Age-dependent population dispersal in RN. Mathematical Biosciences 149, 37–56. International Atomic Energy Agency (IAEA), 2000. Modelling of the Transfer of Radiocaesium from Deposition to Lake Ecosystems. IAEA-TECDOC-1143. International Atomic Energy Agency, Vienna, pp. 111–114. International Atomic Energy Agency (IAEA), 2006. The Chernobyl Forum 2003–2005: Chernobyl’s Legacy: Health Environment and Socio-Economic Impacts and Recommendations to the Governments of Belarus, Second revised version. The Russian Federation and Ukraine. Second revised version, IAEA, Vienna. Joint Norwegian-Russian Group, 1997. Sources Contributing to Radioactive Contamination of the Techa River and Areas Surrounding the “Mayak” Production Association, Urals, Russia. Joint Norwegian-Russian Group for Investigation of Radiation Contamination in the Northern Areas. Norwegian Radiation Protection Authority, Østerås, p. 134. Jørgensen, S.E. (Ed.), 1978. Proceedings of the Conference on Ecological Modelling, International Society for Ecological Modelling. Copenhagen, Denmark, p. 891. Jørgensen, S.E. (Ed.), 1983. Application of Ecological Modelling in Environmental Management, Parts A and B. Elsevier Scientific Publishing Company, Amsterdam-Oxford-New York. Jørgensen, S.E. (Ed.), 1984. Modelling the Fate and Effect of Toxic Substances in the Environment. Elsevier Science Publishers BV, The Netherlands, p. 342. Keats, A., Yee, E., Lien, F.-S., 2007. Efficiently characterizing the origin and decay rate of a nonconservative scalar using probability theory. Ecological Modelling 205, 437–452. Lauenroth, W.K., Skogerboe, G.V., Flug, M. (Eds.), 1983. Analysis of Ecological Systems: State-of-the-Art in Ecological Modelling. Elsevier Scientific Publishing Company, Amsterdam-Oxford-New York, p. 992. Li, B.L., Loehle, C., Malon, D., 1996. Microbial transport through heterogeneous porous media: random walk, fractal and percolation approaches. Ecological Modelling 85, 285–302. Marsh, L.M., Jones, R.E., 1988. The form and consequences of random walk movement models. Journal of Theoretical Biology 133, 113–131. Monte, L., 2002. A methodology for modelling the contamination of moving organisms in water bodies with spatial and time dependent pollution levels. Ecological Modelling 158, 21–33. Monte, L., 2009a. Predicting the effects of ionising radiation on ecosystems by a generic model based on the Lotka-Volterra equations. Journal of Environmental Radioactivity 100, 477–483. Monte, L., 2009b. Multi-model approach and evaluation of the uncertainty of model results, Rationale and applications to predict the behaviour of contaminants in the abiotic components of the fresh water environment. Ecological Modelling 220, 1469–1480. Monte, L., Fratarcangeli, S., Pompei, F., Quagguia, S., Battella, C., 1993. Bioaccumulation of 137 Cs in the main species of fishes in lakes of central Italy. Radiochimica Acta 60, 219–222.

3290

L. Monte / Ecological Modelling 220 (2009) 3280–3290

Odum, E.P., 1959. Fundamentals of Ecology. W.B. Saunders Company, Philadelphia and London, p. 546. O’Toole, Metcalfe, C., Craine, I., Gross, M., 2006. Release of persistent organic contaminants from carcasses of Lake Ontario Chinook salmon (Oncorhynchus tshawytscha). Environmental Pollution 140, 102–113. ˜ Perianez, R., 2007. Chemical and oil spill rapid response modelling in the Strait of Gibraltar-Alborán Sea. Ecological Modelling 207, 210–222. Rajar, R., Cetina, M., Sirca, A., 1997. Hydrodynamic and water quality modelling: case studies. Ecological Modelling 101, 209–228. Rescigno, A., Lambrecht, R., Duncan, C.C., 1983. Mathematical methods in the formulation of pharmacokinetical models. In: Lambrecht, R.M., Rescigno, A. (Eds.), Tracer Kinetics and Physiologic Modeling. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, pp. 59–119. Sazykina, T.G., Kryshev, A.I., 2003a. EPIC database on the effects of chronic radiation in fish: Russian/FSU data. Journal of Environmental Radioactivity 68, 65–87. Sazykina, T.G., Kryshev, I.I., 2003b. Effects of ionising radiation on terrestrial animals: dose–effects relationships. In: International conference on the Protection of the Environment from the Effects of Ionizing Radiation. Contributed papers, Stockholm, Sweden 6–10 October, pp. 95–97. Sazykina, T.G., Kryshev, I.I., 2006. Radiation effects in wild terrestrial vertebrates—the EPIC collection. Journal of Environmental Radioactivity 88, 11–48. Scott, E.M. (Ed.), 2003. Modelling Radioactivity in the Environment. Elsevier Science Ltd, Kidlington, Oxford, UK, p. 427.

Smith, J.T., Comans, R.N.J., 1996. Modelling the diffusive transport and remobilisation of 137 Cs in sediments: the effects of sorption kinetics and reversibility. Geochimica et Cosmochimica Acta 60, 995–1004. Smith, J.T., Wright, S.M., Cross, M.A., Monte, L., Kudelsky, A.V., Saxén, R., Vakulovsky, S.M., Timms, D.N., 2004. Global Analysis of the Riverine Transport of 90Sr and 137 Cs. Environmental Science and Technology 38, 850–857. Turchin, P., 1991. Translating foraging movements in heterogeneous environments into spatial distribution of foragers. Ecology 72, 1253–1266. UNSCEAR, 2000. Report of the United Nations Scientific Committee on the Effects of Atomic Radiation to the General Assembly. Appendix J: Exposures and Effects of the Chernobyl Accident. The United Nation Scientific Committee on the Effects of Atomic Radiation, United Nations. Woodhead, D.S., 2003. A possible approach for the assessment of radiation effects on population of wild organisms in radionuclide-contaminated environment? Journal of Environmental Radioactivity 66, 181–213. Wu, H.i., Li, B.L., Springer, T.A., Neill, W.H., 2000. Modelling animal movement as a persistent random walk in two dimensions: expected magnitude of net displacement. Ecological Modelling 132, 115–134. Zhang, X., Johnson, S.N., Crawford, J.W., Gregory, P.J., Young, I.M., 2007. A general random walk model for the leptokurtic distribution of organism movement: theory and application. Ecological Modelling 200, 79–88. Zheleznyak, M.J., Demchenko, R.I., Khursin, S.L., Kuzmenko, Y.I., Tkalich, P.V., Vitiuk, N.Y., 1992. Mathematical modelling of radionuclide dispersion in the Pripyat–Dnieper aquatic system after the Chernobyl accident. The Science of the Total Environment 112, 89–114.