Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Utilization of bathymetry data to examine lead sediment contamination distributions in Lake Ontario

  • Received: 15 April 2016 Accepted: 20 June 2016 Published: 21 June 2016
  • Bathymetry data offer interesting opportunities for the analysis of contaminant distribution patterns. This research utilized lead surficial sediment sample data from Lake Ontario that were collected by the Canada Centre for Inland Waters in 1968 and 1998. Traditionally, two-dimensional analyses such as dot maps or proportional circle representation have been utilized to examine pollutant levels. Generating area estimates allows for expanded spatial analysis of contaminant distribution patterns. Lake-wide surfaces were derived using the ordinary kriging technique. These were then layered on bathymetry data to examine three-dimensional relationships between observed pollution patterns and lake-bottom features. Spatial variability was observed in both the 1968 and 1998 datasets. Contamination levels in 1998 dropped substantially, especially in areas that were previously the most heavily polluted and above the Probable Effect Level (4660.23 km2 or 26.72% of the common analysis area lake-bottom in 1998 versus 6189.07 km2 or 62.00% in 1968). Conversely, areas below the Threshold Effect Level increased from 922.09 km2 (5.29%) in 1968 to 3484.22 km2 (19.98%) in 1998. In both years, shallow and sill/ridge areas tended to have lower levels of contamination than deeper lake basins or contaminant inflow areas. The 1968 dataset likely provides a more detailed estimation surface as there were more points available for interpolation procedures. The kriging surfaces when combined with bathymetry, sedimentology information, and knowledge of physical processes provide a comprehensive illustration of the contaminant distributions whether they are high (1968) or when loadings are significantly reduced (1998). The results have implications for future sediment assessment programs and survey design on a lake-wide basis. The bathymetry data allowed for enhanced interpretation and an improved understanding of observed lead pollution patterns.

    Citation: K. Wayne Forsythe, Chris H. Marvin, Danielle E. Mitchell, Joseph M. Aversa, Stephen J. Swales, Debbie A. Burniston, James P. Watt, Daniel J. Jakubek, Meghan H. McHenry, Richard R. Shaker. Utilization of bathymetry data to examine lead sediment contamination distributions in Lake Ontario[J]. AIMS Environmental Science, 2016, 3(3): 347-361. doi: 10.3934/environsci.2016.3.347

    Related Papers:

    [1] Yun Kang, Sourav Kumar Sasmal, Komi Messan . A two-patch prey-predator model with predator dispersal driven by the predation strength. Mathematical Biosciences and Engineering, 2017, 14(4): 843-880. doi: 10.3934/mbe.2017046
    [2] Saheb Pal, Nikhil Pal, Sudip Samanta, Joydev Chattopadhyay . Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Mathematical Biosciences and Engineering, 2019, 16(5): 5146-5179. doi: 10.3934/mbe.2019258
    [3] Sangeeta Kumari, Sidharth Menon, Abhirami K . Dynamical system of quokka population depicting Fennecaphobia by Vulpes vulpes. Mathematical Biosciences and Engineering, 2025, 22(6): 1342-1363. doi: 10.3934/mbe.2025050
    [4] Chunmei Zhang, Suli Liu, Jianhua Huang, Weiming Wang . Stability and Hopf bifurcation in an eco-epidemiological system with the cost of anti-predator behaviors. Mathematical Biosciences and Engineering, 2023, 20(5): 8146-8161. doi: 10.3934/mbe.2023354
    [5] Kawkab Al Amri, Qamar J. A Khan, David Greenhalgh . Combined impact of fear and Allee effect in predator-prey interaction models on their growth. Mathematical Biosciences and Engineering, 2024, 21(10): 7211-7252. doi: 10.3934/mbe.2024319
    [6] Yuanfu Shao . Bifurcations of a delayed predator-prey system with fear, refuge for prey and additional food for predator. Mathematical Biosciences and Engineering, 2023, 20(4): 7429-7452. doi: 10.3934/mbe.2023322
    [7] Jinxing Zhao, Yuanfu Shao . Bifurcations of a prey-predator system with fear, refuge and additional food. Mathematical Biosciences and Engineering, 2023, 20(2): 3700-3720. doi: 10.3934/mbe.2023173
    [8] Yuhong Huo, Gourav Mandal, Lakshmi Narayan Guin, Santabrata Chakravarty, Renji Han . Allee effect-driven complexity in a spatiotemporal predator-prey system with fear factor. Mathematical Biosciences and Engineering, 2023, 20(10): 18820-18860. doi: 10.3934/mbe.2023834
    [9] Hongqiuxue Wu, Zhong Li, Mengxin He . Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Mathematical Biosciences and Engineering, 2023, 20(10): 18592-18629. doi: 10.3934/mbe.2023825
    [10] Soumitra Pal, Pankaj Kumar Tiwari, Arvind Kumar Misra, Hao Wang . Fear effect in a three-species food chain model with generalist predator. Mathematical Biosciences and Engineering, 2024, 21(1): 1-33. doi: 10.3934/mbe.2024001
  • Bathymetry data offer interesting opportunities for the analysis of contaminant distribution patterns. This research utilized lead surficial sediment sample data from Lake Ontario that were collected by the Canada Centre for Inland Waters in 1968 and 1998. Traditionally, two-dimensional analyses such as dot maps or proportional circle representation have been utilized to examine pollutant levels. Generating area estimates allows for expanded spatial analysis of contaminant distribution patterns. Lake-wide surfaces were derived using the ordinary kriging technique. These were then layered on bathymetry data to examine three-dimensional relationships between observed pollution patterns and lake-bottom features. Spatial variability was observed in both the 1968 and 1998 datasets. Contamination levels in 1998 dropped substantially, especially in areas that were previously the most heavily polluted and above the Probable Effect Level (4660.23 km2 or 26.72% of the common analysis area lake-bottom in 1998 versus 6189.07 km2 or 62.00% in 1968). Conversely, areas below the Threshold Effect Level increased from 922.09 km2 (5.29%) in 1968 to 3484.22 km2 (19.98%) in 1998. In both years, shallow and sill/ridge areas tended to have lower levels of contamination than deeper lake basins or contaminant inflow areas. The 1968 dataset likely provides a more detailed estimation surface as there were more points available for interpolation procedures. The kriging surfaces when combined with bathymetry, sedimentology information, and knowledge of physical processes provide a comprehensive illustration of the contaminant distributions whether they are high (1968) or when loadings are significantly reduced (1998). The results have implications for future sediment assessment programs and survey design on a lake-wide basis. The bathymetry data allowed for enhanced interpretation and an improved understanding of observed lead pollution patterns.


    For the benefit of both the present and future generations, protecting endangered plants and animals as well as maximizing the preservation of the planet's unique biological resources are the main objectives of biodiversity conservation. An increasing number of empirical findings indicate that human-caused environmental deterioration has impacted natural population habitats [1]. The habitat of populations is fragmented, split into patches and the populations become extinct. Numerous studies [2,3,4] have shown that the creation of ecological corridors between patches is more favorable to the survival of species. Therefore, it is of great importance to study the dispersal behavior of populations.

    In 1977, Freedman and Waltman [5] studied a source–source system with dispersal. Their study indicated that under certain conditions, the total population abundance of dispersed species may be greater than that of non-dispersed species. Holt [6] considered a sink–source system:

    {dN1dt=r1N1+DN2DN1,dN2dt=N2ϕ2(N2)+DN1DN2,

    where Ni represents the population density of the species in the patch i, i=1,2. The first patch is a sink patch which is not self-sustaining and the second patch is a source. r1 represents the species' death rate in patch 1. Here, the dispersal rate D is symmetric. ϕ2(N2) is the per-capita growth rate of the species in patch 2, which monotonically declines due to intraspecific competition. Holt showed that if the death rate in the sink patch and the dispersal rate are small, the species persists in the two patches. Otherwise, the species went to extinction in both patches. After that, more and more scholars studied the different types of diffusion between patches. Wang et al. [7,8,9] did a series of work on two-species competitive systems with one-species' diffusion between patches. In [10], the authors analyzed linear dispersal (the diffusion coefficient is constant) in predator-prey systems with Holling II functional response as follows.

    {dMdt=rM+a12MN1b+N1,dN1dt=r1N1(1N1K1)a21MN1b+N1+D2N2D1N1,dN2dt=r2N2+D1N1D2N2.

    Their study showed that there exists an optimal dispersal that drives the predator into extinction and makes the prey reach the maximal abundance. These results are biologically important when preserving endangered species. On the basis of [10], Xia et al. [11] proposed and studied a predator-prey model with the prey's dispersal and counter-attack behavior as follows.

    dN1dt=N1(r11+kMδγN1)aN1M1+bN1+D2N2D1N1,dN2dt=r2N2+D1N1D2N2,dMdt=saN1M1+bMdMcN1M.

    Here k is the fear factor and c is the counter-attack rate. Compared with a single strategy, the multiple anti-predator strategies are more beneficial to the persistence and the population density of prey. We notice that in [10] and [11], the dispersal of the prey between the two patches is linear. Actually, there are several scholars who investigated the non-constant case, i.e., nonlinear dispersal. For example, Xia et al. [12] proposed a two-patch model with the Allee effect and nonlinear dispersal:

    {dudt=u(ruA+udbu)+Du(vu),dvdt=v(acv)+Dv(uv),

    where A is the Allee effect constant. D is the dispersal coefficient and dispersal is affected by population density. The authors found that the diffusion rate has a certain influence on the solution of the system. When dispersal is high, nonlinear dispersal is more favorable to population survival than linear dispersal. Kang et al. [13] formulated a Rosenzweig-MacArthur prey-predator two-patch model with mobility only in the predator and assumed that the predator moves toward patches with more concentrated prey-predator interactions:

    {dx1dt=r1x1(1x1K1)b1x1y11+b1h1x1,dy1dt=c1b1x1y11+b1h1x1δ1y1+(b1x1y11+b1h1x1ρ21y2b2x2y21+b2h2x2ρ12y1),dx2dt=r2x2(1x2K2)b2x2y21+b2h2x2,dy2dt=c2b2x2y21+b2h2x2δ2y2+(b2x2y21+b2h2x2ρ12y1b1x1y11+b1h1x1ρ21y2),

    where ri is the intrinsic growth of prey at patch i; ki is the prey carrying capacity in patch i; bi is the predator attacking rate in patch i; hi is the predator handling time in patch i; ci is the energy conversion rate in patch i; and δi is the mortality of the predators in patch i(i=1,2). In the presence of dispersal, the dispersal rate of the predators from patch i to patch j depends on the prey-predator interaction term bjxjyj1+bjhjxj in patch j, and the dispersal term of predators from patch i to patch j is modeled by bjxjyj1+bjhjxjρijyi which gives the net dispersal of the predator at patch i as follows: bjxjyj1+bjhjxjρijyibixiyi1+bihixiρjiyj, where ρij is the dispersal constant of predators in patch j moving to patch i. The authors showed that dispersal may stabilize or destabilize the coupled system and dispersal may generate multiple interior equilibria that lead to rich bistable dynamics or may destroy interior equilibria that lead to the extinction of predators in one patch or both patches. In 2019, Mai [14] proposed the following predator-prey metapopulation model over an arbitrary number of patches:

    {dhi(t)dt= εhi(t)(1hi(t)κpi(t)1+hi(t))1jn,jid(αγpi(t)1+hi(t)+(1α))hi(t)+1jn,jid(αγpj(tτ)1+hj(tτ)+(1α))hj(tτ),dpi(t)dt= hi(t)pi(t)1+hi(t)μpi(t),

    where hi and pi denote the densities of prey and predators in patch i (i=1,2,...,n), respectively. n is the total number of patches, d>0 is the rescaled dispersal rate; α[0,1] denotes the fraction of dispersal due to predation avoidance; ε>0 and κ>0 are the rescaled intrinsic growth rate and carrying capacity of the prey, respectively; and a Holling type II functional response function is used where γ>0 is the scaling factor associated with the maximum consumption rate and yield constant. The prey's dispersal is due to random effect and avoidance from predators. The authors' results showed that the delay and the patterns of prey dispersal jointly affect the stability of predator-prey metacommunities and can induce multiple stability switches. For more relevant works, one can see [15,16,17].

    In natural ecosystems, the interaction between predator and prey is one of the most important roles, the predator-prey model has been extensively studied [18,19,20]. Some scholars [21,22] have argued that the mere presence of a predator may alter the prey's behavior and physiology. Its presence has even more impact than direct predation. In 2016, Wang et al. [23] first proposed the following predator–prey model with a fear effect on prey:

    {dudt=r0uf(k,v)duau2puv,dvdt=cpuvmv,

    where f(k,v)=11+kv and the parameter k reflect the degree of the fear effect on the prey from the predators. Furthermore, for the above system, the authors introduced the Holling type II functional response and concluded that the fear effect can lead to population oscillation from supercritical Hopf bifurcation to subcritical Hopf bifurcation. In recent years, more and more works have focused on the fear effect. Liu et al. [24,25] proposed a Leslie–Gower model with the Allee effect on prey and fear effect on the predator. Their study demonatrates that compared with the case without a fear effect on the predator, the new model has richer dynamic behavior and various kinds of bifurcations will occur. About the impact of the fear effect on dynamic behaviors of different ecosystems, one can refer to [26,27,28] for details. Sarkar et al. [29] thought that the fear effect function proposed by Wang et al. [23] was not consistent with the ecological status. The authors introduced the fear effect function f(k,η,P)=η+k(1η)k+P and proposed the following predator–prey model:

    {dNdt=r0N[η+k(1η)k+P]δNγN2βNP1+ξN,dPdt=(θβN1+ξNd)P,

    where k represents the level of the fear and η[0,1] indicates the minimum cost of the fear effect. β, ξ, d, and θ are the per-capita consumption rate of the predator, the processing time to capture each prey, the mortality rate of the predator, and the conversion rate of prey biomass to predator biomass, respectively. The authors obtained that the fear effect can stabilize the predator-prey system at an interior steady state, or it can create the oscillatory coexistence of all the populations. f(k,η,P) is a monotone increasing function with respect to k. Namely, as the level of fear effect k increases, the prey reproduction rate increases. From the ecological point of view, the above fear effect function is not reasonable. Dong et al. [30] proposed a new fear function f(k,η,y)=η+1η1+ky and created a new model as follows.

    {dxdt=r0x(η+1η1+kαy)δxγx2βxy1+ξαx,dydt=(θβx1+ξαxd)y,

    where α>0 is the predator-taxis sensitivity, which has a positive effect on prey and a negative effect on predator population. For the fear function f(k,η,y)=η+1η1+ky, it is easy to obtain that f(k,y)η as k or y. This is why the authors call the parameter η the saturated fear cost. The authors found that the fear effect from predators does not cause the extinction of the prey population. Even if the fear effect is sufficiently large, the prey can survive under the saturated fear cost. Furthermore, when choosing the fear level k as a bifurcation parameter, the model will arise a Hopf bifurcation point. However, when the predator-taxis sensitivity α is chosen as a bifurcation parameter, the model will arise two Hopf bifurcation points. The predator-taxis sensitivity determines the success or failure of the predator invasion under an appropriate fear level. After that, Meng et al. [31] investigated a predator-prey system incorporating asymmetric dispersal and the fear effect, which influences the birth and death rates of the prey species as follows.

    {dx1dt=r1x11+kyd1[1+η(111+ky)]x1γx21αx1y+D2x2D1x1,dx2dt=r2x2+D1x1D2x2,dydt=sαx1yd2y.

    The authors revealed that the fear effect and its maximum cost exhibit significant negative impacts on predator abundance, but they do not influence the density of the prey.

    Through the above literature survey, we know that the population dispersal considered in the patchy model may be related with multiple factors. It is often recognized that when there is a fear effect from the predator, the prey is more likely to scatter. For example, in Yellowstone National Park, due to the fear of being preyed upon by wolves, the elk evades the open river valleys (where wolves are most active in hunting), and migrate to areas such as the forest edges or steep mountains where it is difficult for wolves to track them. Furthermore, the dispersal rate increases with the degree of the fear effect. So what happens to the dynamic behavior when we consider the case where prey dispersal is affected by the fear effect? In this paper, we propose the following predator-prey model with dispersal subject to the fear effect in a two-patch environment:

    {dx1dt=r1x1(1x1K1)a1x1y+D2x2D1[1+η(111+ky)]x1,dx2dt=r2x2+D1[1+η(111+ky)]x1D2x2,dydt=a2a1x1ydy, (1.1)

    where xi represents the density of species x in patch i (i=1,2), and y represents the density of the predator in patch 1. r1 is the intrinsic growth rate of species x in patch 1, and r2 is the death rate of species x in patch 2. K1 is the carrying capacity of species x in patch 1. a1 is the maximal per-capita consumption rate, i.e., the maximum number of prey that can be eaten in each unit time. a2 denotes the conversion of nutrients into the reproduction of the predators. d and e are the natural mortality rate and the intra-species competition coefficient of species y, respectively. D1(D2) is the dispersal rate of species x from patch 1(2) to patch 2(1). It is worth noting that the dispersal rate of the prey increases as predator biomass increases. The dispersal rate deduced by the fear effect from the predator is not likely to increase indefinitely; instead, it may reach a finite value. Thus, we consider the dispersal rate of prey by multiplying a function f(k,η,y)=1+η(111+ky), where k is the level of the fear effect and η is the maximum cost of the fear effect.

    To simplify the parameters of system (1.1), the following dimensionless quantities are applied:

    ˉx1=x1K1, ˉx2=x2K1, ˉy=yK1, ˉt=r1t

    and

    ˉα1=a1K1r1, ˉr2=r2r1, ˉDi=Dir1, ˉk=kK1, ˉα2=a2a1K1r1, ˉd=da2a1K1.

    With bars dropped, system (1.1) can be rewrittten as

    {dx1dt=x1(1x1)α1x1y+D2x2D1[1+η(111+ky)]x1,dx2dt=r2x2+D1[1+η(111+ky)]x1D2x2,dydt=α2y(x1d). (1.2)

    All parameters in system (1.2) are positive and 0<d<1. The initial conditions are

    x1(0)0,x2(0)0,y(0)0. (1.3)

    To the best of our knowledge, this is the first study to discuss a patchy model in which the migration is induced by the fear effect of the predator. Unlike most studies in which the fear effect influences the birth rate of prey, we will investigate how the fear effect from the predator impacts the diffusion rate of prey populations. In fact, our results obtained in this paper can be seen as an important expansion and supplement to those in [13] and [14].

    The rest of the paper is organized as follows. In Section 2, we obtain the permanence of the system (1.2). Existence and local stability of the equilibria are given in Section 3 and Section 4 demonstrates the global stability. In Section 5, we show that system (1.2) undergoes transcritical bifurcation. Also, some numerical simulation results are shown to verify the analytical results in Section 6. Finally, we give a brief discussion in Section 7.

    Theorem 2.1. All solutions of system (1.2) with initial conditions (1.3) are nonnegative and ultimately bounded for t>0.

    Proof. According to Proposition A.17 of [32], it is not difficult to obtain that all the solutions of system (1.2) are always nonnegative. Next, we define a function

    w(t)=x1(t)+x2(t)+α1α2y(t). (2.1)

    Then, by differentiating (2.1) with respect to t, we obtain

    dwdt=x1(1x1)r2x2dα1yx1(2x1)ˉr(x1+x2+α1α2y)1ˉrw,

    where ˉr=min{1,r2,dα2}. Thus, we obtain

    lim supt+w(t)1ˉr:=G0.

    So all solutions of system (1.2) with initial conditions (1.3) are ultimately bounded. This completes the proof.

    Theorem 2.2. System (1.2) is uniformly persistent if 0<D1<D+, where D+=(1dα1G0)(α1+kα2G0)α1+(η+1)kα2G0, α1<1dG0.

    Proof. From (2.1), for ε>0 small enough, there is T>0 such that, for t>T, we have

    x1(t), x2(t)G0+ε:=Gε, y(t)α2α1G0+ε:=Gε.

    According to the first equation of system (1.2),

    dx1dtx1(1α1GεD1[1+η(111+kGε)]x1).

    Setting ε0 in the above inequality leads to

    lim inft+x1(t)1α1G0D1[1+η(111+kα2α1G0)]:=m1.

    Similarly, according to the third equation of system (1.2),

    dydtα2y(m1dy).

    Then we have

    lim inft+y(t)m1d:=m2.

    Finally, it follows from the second equation of system (1.2) and the above results that

    dx2dtD1[1+η(111+km2)]m1(r2+D2)x2.

    We get

    lim inft+x2(t)D1m1r2+D2[1+η(111+km2)]:=m3.

    When 0<D1<D+, we have m1>d and then m2,m3>0. The above implies that system (1.2) is uniformly persistent. This ends the proof.

    In the following we present five conditions, i.e.,

    (H1):0<D1<(1d)(1+D2r2),(H2):D1=(1d)(1+D2r2),(H3):(1d)(1+D2r2)<D1<1+D2r2,(H4):D1=1+D2r2,(H5):D1>1+D2r2.

    One can get the equilibria by setting the right-hand sides of system (1.2) equal to zero. Obviously, the trivial equilibrium E0(0,0,0) always exists. If (H1), (H2), or (H3) holds, system (1.2) also has a predator-extinction equilibrium E1(ˉx1,ˉx2,0), where ˉx1=1D1r2r2+D2 and ˉx2=D1ˉx1r2+D2.

    For the potential coexistence equilibrium E2(x1,x2,y), we have

    x1=d,x2=D1r2+D2[1+η(111+ky)]d,

    and y is the positive root of the the following equation:

    P1y2+P2y+P3=0, (3.1)

    where

    P1=kα1,P2=α2+kηD1r2r2+D2+k(D1r2r2+D2+d1),P3=D1r2r2+D2+d1.

    If D1(1d)(1+D2r2), one has P2>0 and P30. Then Eq. (3.1) has no positive real root. If (H1) holds, it has a unique positive real root y, where

    y=P2+P224P1P32P1.

    Therefore, we can state the following theorem:

    Theorem 3.1. System (1.2) has a unique co-existence equilibrium E2(x1,x2,y) if (H1) holds.

    In the next section, we will discuss the local stability of all equilibria. The Jacobian matrix at the equilibrium E(x1,x2,y) of system (1.2) is calculated as

    JE=(12x1α1yD1[1+η(111+ky)]D2α1x1kηD1x1(1+ky)2D1[1+η(111+ky)]r2D2kηD1x1(1+ky)2α2y0α2(x1d)).

    We then have the following result:

    Theorem 3.2. The trivial equilibrium E0(0,0,0) is

    (i) unstable if (H1), (H2), or (H3) holds;

    (ii) locally asymptotically stable if (H5) holds;

    (iii) a saddle-node if (H4) holds.

    Proof. For the trivial equilibrium E0(0,0,0), the Jacobian matrix is given by

    JE0=(1D1D20D1r2D2000dα2).

    The corresponding characteristic equation is

    (λ+dα2)(λ2(1D1r2D2)λ+D1r2(r2+D2))=0. (3.2)

    We know that Eq. (3.2) has a positive root and two negative roots if (Hi)(i=1,2,3) holds, thus, E0 is unstable. If (H5) holds, one can verify that the three roots of Eq. (3.2) are negative, thus E0 is locally asymptotically stable. Notice that Eq. (3.2) has two negative roots and a unique zero root when (H4) holds. Applying the Taylor expansion of 11+ky at the origin, system (1.2) can be rewritten as

    {dx1dt=a100x1+a010x2+a200x21+a101x1y+O(|x1,y|)3,dx2dt=b100x1+b010x2+b101x1y+O(|x1,y|)3,dydt=c001y+c101x1y, (3.3)

    where

    a100=D2r2, a010=D2, a200=1, a101=ηkD2+ηkr2+α1r2r2, b100=r2+D2r2, b010=r2D2, b101=(r2+D2)ηkr2, c001=dα2, c101=α2.

    Then we make the following transformation:

    (x1x2y)=(r2D2r2+D20110001)(X1X2Y),

    and system (3.3) becomes

    {dX1dt=ˆa200X21+ˆa110X1X2+ˆa101X1Y+ˆa011X2Y+ˆa020X22+O(|X1,X2,Y|)3,dX2dt=ˆb010X2+ˆb101X1Y+ˆb200X21+ˆb110X1X2+ˆb011X2Y+O(|X1,X2,Y|)3,dYdt=ˆc001Y+ˆc101X1Y+ˆc011X2Y, (3.4)

    where

    ˆa200=(r2+D2)r22D2r2+r22+D2, ˆa110=2D2r2D2r2+r22+D2, ˆa101=(r2+D2)r2(ηk+α1)D2r2+r22+D2,ˆa011=D2(ηk+α1)D2r2+r22+D2, ˆa020=D22(D2r2+r22+D2)(r2+D2), ˆb010=D2r2+r22+D2r2,ˆb101=(ηkr22+(D2ηk+ηk+α1)r2+D2ηk)(r2+D2)D2r2+r22+D2, ˆb200=D22(D2r2+r22+D2)(r2+D2),ˆb110=2D2r2D2r2+r22+D2, ˆb011=(ηkr22+(D2ηk+ηk+α1)r2+D2ηk)D2r2(D2r2+r22+D2),ˆc001=dα2, ˆc101=r2α2, ˆc011=D2α2r2+D2.

    For X10, there exists a center manifold

    X2=ˆb200ˆb010X21+O(|X1|)3, Y=O(|X1|)3.

    Then on the center manifold, system (3.4) becomes

    dX1dt=ˆa200X21+O(|X1|)3.

    Obviously, the coefficient of X21 is ˆa200=(r2+D2)r22D2r2+r22+D2<0. Therefore, E0 is a saddle-node. This ends the proof of Theorem 3.2.

    Theorem 3.3. The predator-extinction equilibrium E1(ˉx1,ˉx2,0) is

    (i) unstable if (H1) holds;

    (ii) locally asymptotically stable if (H3) holds;

    (iii) a saddle-node if (H2) holds.

    Proof. For the predator-extinction equilibrium E1(ˉx1,ˉx2,0), the Jacobian matrix is given by

    JE1=(D1r2D2r2+D21D2(α1+kηD1)(1D1r2r2+D2)D1r2D2kηD1(1D1r2r2+D2)00α2(1D1r2r2+D2d)).

    JE1 always has an eigenvalue λ=α2(1D1r2r2+D2d).

    (i) If (H1) holds, one can verify that λ>0, thus E1(ˉx1,ˉx2,0) is unstable.

    (ii) If (H3) holds, then we have λ<0. Define

    J1=(D1r2D2r2+D21D2D1r2D2).

    The determinant and trace of J1 can be simplified as

    Det(J1)=D1r2+r2+D2,Tr(J1)=D1r2D2r2+D21r2D2.

    Combining with the conditions for the existence of E1, we obtain Det(J1)>0 and Tr(J1)<0. It follows that J1 has two negative eigenvalues. Hence, if (H3) holds, JE1 has three negative eigenvalues which implies that E1 is locally asymptotically stable.

    (iii) If (H2) holds, we also can deduce that Det(J1)>0 and Tr(J1)<0. Now, JE1 has a unique zero eigenvalue λ=0. Using the same method of proof as in Theorem 3.2, we can get that E1 is a saddle-node. This completes the proof.

    Theorem 3.4. When (H1) holds, the co-existence equilibrium E2(x1,x2,y) exists and is always locally asymptotically stable.

    Proof. The Jacobian matrix at E2(x1,x2,y) can be calculated as

    JE2=(12x1α1yD1[1+η(111+ky)]D2α1x1kηD1x1(1+ky)2D1[1+η(111+ky)]r2D2kηD1x1(1+ky)2α2y0α2(x1d)).

    The corresponding characteristic equation is

    λ3+n1λ2+n2λ+n3=0,

    where

    n1= r2+D2+x1+D1D2r2+D2[1+η(111+ky)]>0,n2= α2y[α1x1+kηD1x1(1+ky)2]+(r2+D2)x1>0,n3= (r2+D2)α1α2x1y+kηD1α2r2x1y(1+ky)2>0.

    After a simple calculation, we obtain that

    n1n2n3={x1+D1D2r2+D2[1+η(111+ky)]}n2+r2N+D2(α2ykηD1x1(1+ky)2+(r2+D2)x1)>0.

    Applying the Routh-Hurwitz criteria[33], when E2 exists, i.e., (H1) holds, we can verify that E2 is locally asymptotically stable. This completes the proof.

    The existence and stability conditions for all equilibria are given in Table 1.

    Table 1.  Existence and local stability of all equilibria.
    Equilibrium Existence Stability
    E0(0,0,0) Always (H1), (H2), or (H3), Unstable (H4), Saddle-node (H5), Stable
    E1(ˉx1,ˉx2,0) (H1), (H2), or (H3) (H3), Stable (H2), Saddle-node (H1), Unstable
    E2(x1,x2,y) (H1) Stable

     | Show Table
    DownLoad: CSV

    In this section we shall analyze the global stability of all the equilibria under some conditions. Consider the following equations:

    {dx1dt=x1(1x1)α1x1y+D2x2D1[1+η(111+ky)]x1,dx2dt=r2x2+D1[1+η(111+ky)]x1D2x2. (4.1)

    Let U1=x1 and U2=x1+x2. Then model (4.1) can be transformed into the following system:

    {dU1dt=U1(1U1)(D1+D2)U1+D2U2α1U1yD1η(111+ky)U1,dU2dt=U1(1U1)+r2U1r2U2α1U1y. (4.2)

    By the nonnegativity of the solutions, from Eq. (4.2), we obtain

    {dU1dtU1(1U1)(D1+D2)U1+D2U2,dU2dtU1(1U1)+r2U1r2U2.

    Applying the comparison theorem for differential equations, we establish the comparison equation:

    {dN1dt=N1(1N1)(D1+D2)N1+D2N2,dN2dt=N1(1N1)+r2N1r2N2. (4.3)

    Then, for model (4.3), we can obtain the following lemmas:

    Lemma 4.1. System (4.3) has a trivial equilibrium P0(0,0). If (H1), (H2), or (H3) holds, system (4.3) has a unique positive equilibrium P1(N+1,N+2), where N+1=1D1r2r2+D2,N+2=D1N+1r2+D2+N+1.

    Lemma 4.2. The trivial equilibrium P0(0,0) is globally asymptotically stable if (H5) holds.

    Proof. The Jacobian matrix of model (4.3) at P0(0,0) is given by

    JP0=(1(D1+D2)D21+r2r2).

    The two eigenvalues of JP0 satisfy the following equations:

    λ1+λ2=1r2D1D2,λ1λ2=D1r2(r2+D2).

    If (H5) holds, JP0 has two negative eigenvalues, thus, P0(0,0) is locally asymptotically stable. Furthermore, there is no positive equilibrium when (H5) holds, and then there exists no limit cycle in the first quadrant. So P0(0,0) is globally asymptotically stable.

    Lemma 4.3. The positive equilibrium P1(N+1,N+2) is globally asymptotically stable if (H1), (H2), or (H3) holds.

    Proof. The Jacobian matrix at P1(N+1,N+2) is

    JP1=(D1r2D2r2+D21D2D22D1r2r2+D21+r2r2).

    The two eigenvalues of JP1 satisfy the following equations:

    λ1+λ2=D1r2D2r2+D21r2D2,λ1λ2=(r2+D2)D1r2.

    Therefore, by some simple calculations, if (H1), (H2), or (H3) holds, the two eigenvalues λ1,2<0 hold, and obviously, P1(N+1,N+2) is locally asymptotically stable. Next, from (4.3), we have

    (r2+D2)(N1N2)+d(N1N2)dt=D1N10.

    Hence, as t,

    lim supt+(N1(t)N2(t))0.

    Then, we denote the right-hand side of system (4.3) as F(N1,N2) and G(N1,N2), respectively. Introduce a Dulac function B(N1,N2)=1N1N2. It follows that

    (BF)N1+(BG)N2=D2N211+r2+N2N1N22<0.

    Applying the Bendixson-Dulac discriminant, system (4.3) has no limit cycle in the first quadrant. Therefore, when P1(N+1,N+2) exists, i.e., (H1), (H2), or (H3) holds, P1(N+1,N+2) is globally asymptotically stable.

    According to Lemma 4.2, we can get

    limt+x1(t)=0, limt+x2(t)=0.

    The third equation of system (1.2) can be rewritten as the following limit equation:

    dydt=dα2y.

    Then we can obtain

    limt+y(t)=0. (4.4)

    Combining the above results, we have the following theorem:

    Theorem 4.1. For system (1.2), the trivial equilibrium E0 is globally asymptotically stable if (H5) holds.

    Further, Lemma 4.3 yields that

    lim supt+U1(t)limt+N1(t)=N+1. (4.5)

    Assume that d>1D1r2r2+D2, i.e., (H3), (H4), or (H5) holds. Thus, there is ε0>0 small enough such that 1D1r2r2+D2+ε0d<0. From Eq. (4.5), for ε0>0, there is T>0 such that, for t>T, we have x1(t)<1D1r2r2+D2+ε0. Then the third equation of system (1.2) can be rewritten as

    dydt<α2y(1D1r2r2+D2+ε0d)<0.

    It follows that

    limt+y(t)=0. (4.6)

    Thus, when y(t)0, the limit equation of system (1.2) is

    {dx1dt=x1(1x1)+D2x2D1x1,dx2dt=r2x2+D1x1D2x2.

    Combining with Theorem 2.2 in [9], we obtain

    limt+x1(t)=1r1D2r2+D2=ˉx1, limt+x2(t)=D1ˉx1r2+D2=ˉx2.

    In summary, we can state the following theorem:

    Theorem 4.2. For system (1.2), the predator-extinction equilibrium E1 is globally asymptotically stable if (H3) holds.

    Next, by using the geometric approach [34,35], we will discuss the global stability of the positive equilibrium E2. We denote the variational matrix of the system (1.2) by V as

    V=(12x1α1yD1[1+η(111+ky)]D2α1x1kηD1x1(1+ky)2D1[1+η(111+ky)]r2D2kηD1x1(1+ky)2α2y0α2(x1d)).

    Then the second additive compound matrix of V is given by

    V|2|=(MkηD1x1(1+ky)2α1x1+kηD1x1(1+ky)20ND2α2yD1[1+η(111+ky)]r2D2+α2(x1d)),

    where M=12x1α1yD1[1+η(111+ky)]r2D2 and N=12x1α1yD1[1+η(111+ky)]+α2(x1d). We define a diagonal matrix P = diag{1,x1y,x1y}. It can be obtained that P1 = diag{1,yx1,yx1} and

    Pf=dPdt=diag{0,˙x1yx1˙yy2,˙x1yx1˙yy2}.

    Hence

    PfP1=diag{0,˙x1x1˙yy,˙x1x1˙yy}, PV|2|P1=V|2|.

    Then we state

    A=PfP1+PV|2|P1=(a11a12a21a22),

    where

    a11=12x1α1yD1[1+η(111+ky)]r2D2,a12=(kηD1x1(1+ky)2,α1x1+kηD1x1(1+ky)2),a21=(0,α2y)T,a22=(BD2D1[1+η(111+ky)]C),B=˙x1x1˙yy+12x1α1yD1[1+η(111+ky)]+α2(x1d),C=˙x1x1˙yyr2D2+α2(x1d).

    Next, let (u,v,w) be a vector in R3 and its norm |.| is defined as |(u,v,w)|=max{|u|,|v+w|}. Denote the Lozinski measures with respect to this norm and with respect to the L1 norm by μ and μ1, respectively. Then

    μ(A)sup(g1,g2)sup{μ1(a11)+|a12|,μ1(a22)+|a21|}

    where |a12|,|a21| are matrix norms with respect to the L1 vector norm. We have

    μ1(a11)= 12x1α1yD1[1+η(111+ky)]r2D2,|a12|= α1x1+kηD1x1(1+ky)2, |a21|=α2y,μ1(a22)= ˙x1x1+max{12x1α1y,r2}.

    Using equations of system (1.2) and with further simplification, we get the following:

    g1=μ1(a11)+|a12|=˙x1x1(1α1kηD1(1+ky)2)x1r2D2D2x2x1˙x1x1(r2+D2)(1α1kηD1(1+ky)2)x1,g2=μ1(a22)+|a21|=˙x1x1+α2y+max{12x1α1y,r2}=˙x1x1+max{12x1(α1α2)y,r2+α2y}.

    Therefore, we have

    μ(A)˙x1x1min{m,n},

    where

    m= r2+D2+(1α1kηD1(1+ky)2)x1,n= max{12x1(α1α2)y,r2+α2y}.

    Now, we denote γ1=lim inft+{x1(t),x2(t),y(t)} and γ2=max{1(2+α1α2)γ1,r2+G0α22α1}, where G0=lim supt+y(t). So we can deduce that

    μ(A)˙x1x1min{δ1,γ2},

    where

    δ1=r2+D2+(1α1kηD1(1+kγ1)2)γ1.

    Then we define

    γ3=min{δ1,γ2}.

    Thus, by the above discussion, we state

    μ(A)˙x1x1γ31tt0μ(A)ds1tln|x2(t)x2(0)|γ3.

    Therefore, if γ3>0, we have

    lim supt+1tt0μ(A)dsγ3<0.

    Then we can state the following theorem:

    Theorem 4.3. For system (1.2), the co-existence equilibrium E2 is globally asymptotically stable if (H1) holds and γ3>0, where γ3=min{δ1,γ2}, γ2=max{1(2+α1α2)γ1,r2+G0α22α1} and γ1=lim inft+{x1(t),x2(t),y(t)}.

    Remark 4.1. Specially, when D11α1kη, α1<1, γ1>12+α1α2, α22<α1r2G0, the conditions of Theorem 4.3 hold. In other words, when the prey's dispersal from the source to the sink (i.e.,D1) is small, the maximal per-capita consumption rate and the conversion of nutrients into the reproduction of the predators are not large, all species in both patches can coexist. Through numerical simulations in Section 6, we find that when E2 exists, i.e., (H1) holds, the uinque positive equilibrium is globally asymptotically stable. However, we cannot present the rigorous proof. We will investigate the above in future work.

    From Section 3, we can see that the predator-extinction equilibrium E1(1D1r2r2+D2,D1(1D1r2r2+D2)r2+D2,0) coincides with the trivial equilibrium E0(0,0,0) if (H4) holds. In the following, we will investigate the existence of transcritical bifurcation in system (1.2) by using the standard Sotomayor theorem [36].

    Theorem 5.1. System (1.2) experiences a transcritical bifurcation around the trivial equilibrium E0 when D1=D1=1+D2r2 and D2r221r2.

    Proof. Notice that the Jacobi matrix JE1 possesses a unique zero eigenvalue when D1=D1. Let

    α:=(α1α2α3)=(r210), β:=(β1β2β3)=(r2+D2D210)

    be the eigenvectors of JE1 and JTE1 corresponding to a zero eigenvalue, respectively. Denote

    F=(F1F2F3)=(x1(1x1)α1x1y+D2x2D1[1+η(111+ky)]x1r2x2+D1[1+η(111+ky)]x1D2x2α2y(x1d)).

    Then we have

    FD1(E0;D1)=(000),
    DFD1(E0;D1)α=(r210),
    D2F(E0;D1)(α,α)=(2F1x21α21+22F1x1x2α1α2+2F1x22α22+22F1x2yα2α3+2F1y2α23+22F1x1yα1α32F2x21α21+22F2x1x2α1α2+2F2x22α22+22F2x2yα2α3+2F2y2α23+22F2x1yα1α32F3x21α21+22F3x1x2α1α2+2F3x22α22+22F3x2yα2α3+2F3y2α23+22F3x1yα1α3)(E0;D1)=(2r2200).

    It follows from above that

    βTFD1(E0;D1)=0,βTDFD1(E0;D1)α=1r2(r2+D2)D2,βTD2F(E0;D1)(α,α)=2r22(r2+D2)D20.

    Therefore, according to the Sotomayor theorem, we know that system (1.2) experiences a transcritical bifurcation at E0 if D2r221r2. Theorem 5.1 is proved.

    To validate the theoretical findings in the preceding parts, we will present the MATLAB numerical simulations of system (1.2) in this section.

    Case I. Stability of the equilibria

    Here we take the parameters of system (1.2) as follows:

    α1=0.4, η=0.2, k=0.15, α2=0.3, d=0.5. (6.1)

    In Figures 13, the equilibria E0(0,0,0), E1(0.25,0.625,0), and E2(0.5,0.6319969406,0.3952971308) are globally asymptotically stable, respectively. The ecological significance of these results are well understood as follows. First, when the prey's dispersal from the source to the sink is large (i.e.,D1>1+D2r2), all populations will be extinct (see Figure 1). Second, when the dispersal is intermediate, i.e., (1d)(1+D2r2)<D1<1+D2r2, the predator goes to extinction and the prey persists in the two patches (see Figure 2). The total prey population abundance is T(D1)=(1D1r2r2+D2)(1+D1r2+D2). Figure 4 shows the trajectory of T about D1, where D1=1r22r2,D01=2D1, and Tmax=T(D1). If 0<D1<D01, the population density of the prey with diffusion is higher than that without diffusion. Furthermore, we can infer that there is an optimal diffusion coefficient D1 that maximizes the total population density of the prey. Third, when the dispersal is small (i.e.,0<D1<(1d)(1+D2r2)), the species in both patches will be permanent (see Figure 3).

    Figure 1.  Global stability of E0 for D1=1.5, D2=0.1, and r2=0.3 with other parameters chosen in (6.1).
    Figure 2.  Global stability of E1 for D1=1, D2=0.1, and r2=0.3 with other parameters chosen in (6.1).
    Figure 3.  Global stability of E2 for D1=0.5, D2=0.1, and r2=0.3 with other parameters chosen in (6.1).
    Figure 4.  The curve of the total prey population density T about D1.

    Case II. Impact of the fear effect on system (1.2)

    For this purpose, we take the following parameters:

    α1=0.4, r2=0.3, α2=0.3, d=0.5, η=0.2, k=0.15, D1=0.5, D2=0.1. (6.2)

    When the dispersal is small (e.g.,0<D1<(1d)(1+D2r2)), the prey and predator coexist at the equilibrium E2(x1,x2,y). Then we know that the fear effect k has no impact on the existence and local stability of the co-existence equilibrium. It only influences the population densities of species x2,y. In Figure 5, we present the dynamic behavior of x1, x2, and y for different levels of fear effect k over time t, respectively. We observe that as k increases, the steady state level of y gradually decreases, but the steady state level of x2 will increase. In Figure 6, we give the dynamic behavior of the populations for different values of η, i.e., the maximum cost of fear. Similarly, we find that as η increases, the steady state level of y will decrease and the steady state level of x2 will increase. In other words, from Figures 56, we obtain that the fear effect plays a negative role for specie y in the first patch. The higher the cost of the prey's fear effect from the predator, the lower the density will be of the population y. In contrast, for specie x in the two patches, the fear effect has a beneficial impact on the total density.

    Figure 5.  Trajectories of x1, x2, and y with respect to time t for different values of k and other parameters chosen in (6.2).
    Figure 6.  Trajectories of x1, x2, and y with respect to time t for different values of η and other parameters chosen in (6.2).

    In this paper, we have analyzed a patchy model in which the migration is induced by the fear effect from the predator. We have investigated that prey species in the source patch can migrate to the sink patch at a faster rate due to the fear effect. Meanwhile, the dispersal rate induced by the fear effect does not increase indefinitely; instead, it may reach a finite value. In other words, we have introduced a fear effect factor k and a maximum fear cost η in our model.

    We have explored the nonnegativity and boundedness of model (1.2). Also, we have obtained that the model is uniformly persistent when 0<D1<D+. For system (1.2), it always has a trivial equilibrium. There exists a predator-extinction equilibrium if (H1), (H2), or (H3) holds and a unique co-existence equilibrium if (H1) holds. By calculating the eigenvalues of the Jacobian matrices at each equilibrium, and combining with the central mainfold theorem, we have obtained the local stability of all equilibria. Then, using the comparison theorem and Dulac's theorem, we have obtained the global stability of the boundary equilibrium. Moreover, a geometric approach in [34,35] has been used to prove the global stability of a unique co-existence equilibrium E2(x1,x2,y). Choosing the diffusion coefficient D1 as the bifurcation parameter, we can conclude that system (1.2) undergoes a transcritical bifurcation. The biological implication of the existence of transcritical bifurcation is that, when (H1), (H2), or (H3) holds, the predator goes extinct, but the species x still exists. If D1 increases, species x easily migrates to patch 2 and will go extinct.

    After that, we have verified the obtained theoretical results by performing numerical simulations through MATLAB. We have selected suitable parameters and analyzed the global stability of all equilibria, which is in agreement with our theoretical results. When the dispersal is intermediate, we find an optimal diffusion coefficient to make the density of the prey population reach its maximum. In Figures 56, we have explored the impacts of fear effect k and maximum fear cost η on the dynamic behavior of system (1.2). They play a negative role for specie y in the first patch and have a beneficial impact on the total density of the prey.

    Unlike [13], we consider the prey's dispersal instead of the predator's disperal between patches. Also, unlike [14], the prey's dispersal rate is affected by the fear effect from the predator in our paper. This is the first time that this type of dispersal has been considered. Overall, this paper extends the literatures [13,14] by considering the impact of dispersal induced by the fear effect. In our model (1.2), dispersal plays a significant role in a the persistence of populations in a source-sink patch. Furthermore, our results show that the fear effect is beneficial to the total population density of the prey, which is different from the results that the fear effect has no influence on the density of prey in [31]. This is because we posit that fear does not directly influence the birth and death rates of the prey species but instead induces prey migration as a response to predation risk. In other words, the prey successfully evades the predator, thereby the dispersal deduced by the fear effect contributes to an increase in prey population density. Our research has important ramifications for protecting threatened and endangered species. By creating ecological corridors and regulating the rate of species movement across patches, we can lessen the consequences of isolated patches and eventually raise the overall abundance of endangered species.

    Our research shows that the fear effect only impacts the total population density, and it has no influence on the existence and stability of the equilibria. This is possibly explained by the fact that we considered the Holling type I functional response between the prey and the predator. Therefore, in future research, we will consider other functional responses. Furthermore, we have noticed that some scholars focused on inverse problems [37,38,39] and boundary observation [40,41,42] in studying multi-population aggregations. Their goal was to reconstruct the diffusion coefficient, advection coefficient, and interaction kernels of the aggregation system, which characterize the dynamics of different populations. In future work, we will connect our study on inverse problems in mathematical biology with more application backgrounds.

    Jin Zhong: writing-original draft; Lijuan Chen and Fengde Chen: supervision, writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported by the National Natural Science Foundation of China (11601085) and the Natural Science Foundation of Fujian Province (2021J01614, 2021J01613).

    The authors declare there is no conflict of interest.

    [1] U.S. Environmental Protection Agency and Government of Canada (USEPAGC), (1995) The Great Lakes—An Environmental Atlas and Resource Book (3rd Edition), Chicago: Great Lakes National Program Office.
    [2] Great Lakes Information Network (GLIN), 2004. Great Lakes Facts and Figures. Great Lakes Information Network. Available from: http://www.great-lakes.net/lakes/ref/lakefact.html
    [3] Crane J, MacDonald DD (2003) Applications of Numerical Sediment Quality Targets for Assessing Sediment Quality Conditions in a US Great Lakes Area of Concern. Environ Manage 32: 128-140.
    [4] Canadian Council of Ministers of the Environment (CCME), 1999. Canadian Environmental Quality Guidelines. Winnipeg, MB: Canadian Council of Ministers of the Environment.
    [5] Forsythe KW, Marvin CH, Valancius CJ, et al. (2016) Using Geovisualization to Assess Lead Sediment Contamination in Lake St. Clair. The Canadian Geographer/Le Géographe canadien 60: 149-158.
    [6] Newell RG, Rogers K, The U.S. Experience with the phasedown of lead in gasoline, 2003. Available from: http://web.mit.edu/ckolstad/www/Newell.pdf
    [7] Currie E (1994) Contaminated Sediment Removal Program Environment Canada—Great Lakes Cleanup Fund.
    [8] Forsythe KW, Watt JP (2006) Using Geovisualization to Assess Sediment Contamination in the “Sixth” Great Lake. In: Strobl J, Blaschke T, Griesebner G (Eds), Proceedings of the 18th Symposium for Applied Geographic Information Processing, AGIT XVIII, Salzburg, Austria, July 5-7, 2006, Herbert Wichmann Verlag, Hüthig GmbH & Co. KG: Heidelberg, Germany, 161-170.
    [9] Forsythe KW, Marvin CH, Valancius CJ, et al. (2016) Geovisualization of Mercury Contamination in Lake St. Clair Sediments. J Mar Sci Eng 4: 19.
    [10] Forsythe KW, Dennis M, Marvin CH (2004) Comparison of Mercury and Lead Sediment Concentrations in Lake Ontario (1968–1998) and Lake Erie (1971-1997/98) using a GIS-based Kriging Approach. Water Qual Res J Can 39: 190-206.
    [11] Jakubek DJ, Forsythe KW (2004) A GIS-based kriging approach for assessing Lake Ontario sediment contamination. Great Lakes Geographer 11: 1-14.
    [12] Forsythe KW, Marvin CH (2005) Analyzing the Spatial Distribution of Sediment Contamination in the Lower Great Lakes. Water Qual Res J Can 40: 389-401.
    [13] Gewurtz SB, Helm PA, Waltho J, et al. (2007) Spatial distributions and temporal trends in sediment contamination in Lake St. Clair. J Great Lakes Res 33: 668-685.
    [14] Gewurtz SB, Shen L, Helm PA, et al. (2008) Spatial distribution of legacy contaminants in sediments of Lakes Huron and Superior. J Great Lakes Res 34: 153-168.
    [15] Forsythe KW, Marvin CH (2009) Assessing Historical versus Contemporary Mercury and Lead Contamination in Lake Huron Sediments. Aquat Ecosyst Health 12: 101-109.
    [16] Forsythe KW, Paudel K, Marvin CH (2010) Geospatial analysis of zinc contamination in Lake Ontario sediments. J Environ Inform 16: 1-10. doi: 10.3808/jei.201000172
    [17] Gawedzki A, Forsythe KW (2012) Assessing anthracene and arsenic contamination within Buffalo River sediments. Int J Ecol 2012: Article ID 496740.
    [18] Forsythe KW, Irvine KN, Atkinson DM, et al. (2015) Assessing lead contamination in Buffalo River sediments. J Environ Inform 26: 106-111.
    [19] Johnston K, Ver Hoef J, Krivoruchko K, et al. (2001) Using ArcGIS geostatistical analyst. Redlands, CA: Environmental Systems Research Institute.
    [20] Isaaks EH, Srivastava MR (1989) An Introduction to Applied Geostatistics. New York, NY: Oxford University Press.
    [21] National Oceanographic and Atmospheric Administration (NOAA), Great Lakes bathymetry. 2014. Available from: http://www.ngdc.noaa.gov/mgg/greatlakes/greatlakes.html.
    [22] Great Lakes Information Network (GLIN), 2004. Lake Ontario Facts and Figures. Great Lakes Information Network. Available from: http://www.great-lakes.net/lakes/ref/ontfact.html
    [23] Marvin CH, Charlton MN, Reiner EJ, et al. (2002) Surficial sediment contamination in Lakes Erie and Ontario: A comparative analysis. J Great Lakes Res 28: 437-450. doi: 10.1016/S0380-1330(02)70596-0
    [24] Marvin CH, Painter S, Rossmann R (2004) Spatial and temporal patterns in mercury contamination in sediments of the Laurentian Great Lakes. Environ Res 95: 351-362.
    [25] Marvin CH, Painter S, Charlton MN, et al. (2004) Trends in spatial and temporal levels of persistent organic pollutants in Lake Erie sediments. Chemosphere 54: 33-40. doi: 10.1016/S0045-6535(03)00660-X
    [26] Dunn RJK, Zigic S, Burling M, et al. (2015) Hydrodynamic and Sediment Modelling within a Macro Tidal Estuary: Port Curtis Estuary, Australia. J Mar Sci Eng 3: 720-744. doi: 10.3390/jmse3030720
    [27] Hessl A, Miller J, Kernan J, et al. (2007) Mapping Paleo-Fire Boundaries from Binary Point Data: comparing Interpolation Methods. The Professional Geographer 59: 87-104. doi: 10.1111/j.1467-9272.2007.00593.x
    [28] Niu Y, Jiao W, Yu H, et al. (2015) Spatial Evaluation of Heavy Metals Concentrations in the Surface Sediment of Taihu Lake. Int J Environ Res Public Health 12: 15028-15039. doi: 10.3390/ijerph121214966
    [29] Ouyang Y, Higman J, Campbell D, et al. (2003) Three-Dimensional Kriging Analysis of Sediment Mercury Distribution: A Case Study. J Am Water Resour Assoc 39: 689-702. doi: 10.1111/j.1752-1688.2003.tb03685.x
    [30] Marvin CH, Charlton MN, Stern GA, et al. (2003) Spatial and temporal trends in sediment contamination in Lake Ontario. J Great Lakes Res 29: 317-331. doi: 10.1016/S0380-1330(03)70437-7
    [31] Beletsky D, Saylor JH, Schwab DJ (1999). Mean circulation in the Great Lakes. J Great Lakes Res 25: 78-93. doi: 10.1016/S0380-1330(99)70718-5
    [32] Environment Canada, Environment Canada’s gasoline regulations, 2009. Available from: http://www.ec.gc.ca/lcpe-cepa/default.asp?lang¼En&n¼54FE5535-1&wsdoc¼8E3C2E9B-38A8-461A-8EC3-C3AA3B1FD585.
    [33] Marvin CH, Charlton MN, Stern GA, et al. (2003) Spatial and temporal trends in sediment contamination in Lake Ontario. J Great Lakes Res 29: 317-331. doi: 10.1016/S0380-1330(03)70437-7
  • Reader Comments
  • © 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(7105) PDF downloads(1254) Cited by(4)

Figures and Tables

Figures(10)  /  Tables(6)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog