A finite element method to solve a population dynamics stage-structured model of intertidal barnacles

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

available at www.sciencedirect.com

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

A finite element method to solve a population dynamics stage-structured model of intertidal barnacles Ana Paula C. Rio Doce a,∗ , Regina C. Almeida b,1 , Michel I.da S. Costa c,2 a

´ ˜ Cient´ıfica-LNCC/MCT, Brazil The Computational Modelling Program, Laboratorio Nacional de Computac¸ao ´ ˜ Cient´ıfica-LNCC/MCT, Av. Getulio ´ Vargas, Department of Computational Mechanics, Laboratorio Nacional de Computac¸ao ´ 333-Petropolis, RJ CEP 25651-075, Brazil c Department of Systems and Control, Laboratorio ´ ˜ Cient´ıfica-LNCC/MCT, Av. Getulio ´ Vargas, 333-Petropolis, ´ Nacional de Computac¸ao RJ CEP 25651-075, Brazil b

a r t i c l e

i n f o

a b s t r a c t

Keywords:

Many marine organisms are fixed or highly sedentary as adults but the adult population

Population dynamics

may be strongly dependent on the oceanic transport of planktonic larvae. In order to assess

Finite element methods

interactions between oceanographic and biological processes that determine the population

Stage-structured life cycle

dynamics of marine organisms with a sessile adult phase restricted to the coastline and a

Biological and physical coupling

planktonic larval phase, we present a stage-structured finite element model for the barnacle Balanus glandula that inhabits the rocky intertidal zone of central California, USA. The larval dispersal depends on knowing the flow pattern. To this end, we develop a numerical procedure to couple physical and biological models in a very simple way when the velocity flow field is known on some discrete values in the domain of interest and at a given time. We investigate the effects of different flow patterns and velocity speeds on the abundance and the distribution of this organism. © 2008 Elsevier B.V. All rights reserved.

1.

Introduction

Most coastal benthic invertebrates have a complex life cycle in which the adult phase is preceded by a planktonic larval phase. For such species, recruitment is one of the main processes governing the temporal fluctuations and the spatial structure of adult population. The recruitment depends on a variety of physical and biological factors, which include spawning, larval dispersal and survival, larval settlement, metamorphosis, post-settlement events, interspecific and intraspecific competitions. In particular, when mature organisms are fixed or highly sedentary, the dispersal in the early life stages is a critical aspect of their population dynamics

∗

(Largier, 2003). During these stages, the larvae are typically planktonic with limited or no swimming abilities so as they are passively transported by the currents. This indicates that larval dispersal depends on coastal oceanographic processes. Therefore, ocean flows affect abundance and distribution patterns through their influence on recruitment processes (Gaylord and Gaines, 2000). Researches in this area have important implications for understanding population dynamics, species geographical ranges and management, spread of invading species and the design of marine reserves (Pineda, 2000). The biological value of recognizing the correlation between near-shore transport in the water and coastal recruitment motivates the development of models for the population dynamics which combine biologi-

Corresponding author. Tel.: +55 24 2233 6135; fax: +55 24 2233 6165. E-mail addresses: [email protected] (A.P.C. Rio Doce), [email protected] (R.C. Almeida), [email protected] (M.I.daS. Costa). 1 Tel.: +55 24 2233 6155; fax: +55 24 2233 6165. 2 Tel.: +55 24 2233 6078; fax: +55 24 2233 6165. 0304-3800/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2008.01.008

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

cal and physical processes (Alexander and Roughgarden, 1996; Gaylord and Gaines, 2000; Ramzi et al., 2006). Roughgarden et al. (1985) introduced a demographic model for a local population of sessile marine invertebrates that have a pelagic larval phase. The processes considered in this model were the settling of larvae onto empty space, and growth and mortality of the settled organisms. The model predicts the number of animals of each stage in the local system through time. Later, Roughgarden et al. (1988) presented a stage model that included mesoscale processes related to the transport of larvae. The other processes included in the model were settlement, mortality and reproduction of adult, and mortality and dispersion of larvae. In this version, larvae are produced by the adults, transported in the water column by advection (only alongshore flow) and eddy-diffusion, and may settle where they strike the coastline. On settling, it is assumed that larvae metamorphose into adults immediately. The model shows that the adult population initially distributed on a small area can spread to adjacent suitable habitats, through the larval dispersion. Possingham and Roughgarden (1990) more fully developed the previous work, showing how the velocity speed of the alongshore flow field, the amount of available suitable habitat and the initial conditions interact to influence species persistence. Later, Alexander and Roughgarden (1996) introduced the hypothesis that recruitment pulses result from the approach and eventual collision of upwelling fronts with the intertidal zone. The results showed how the interaction of offshore and alongshore advection, eddy-diffusion, and the location of the offshore front determine the dynamics of coastal barnacle populations (Alexander and Roughgarden, 1996). In 2000, Gaylord and Gaines extended the model in Roughgarden et al. (1988) so as to better explore the ocean circulation pattern in the geographical distribution of the marine species that have a planktonic life-history phase. They studied four pre-set current patterns: alongshore flow, converging flow, diverging flow and eddy-circulation. Moreover, they considered that the larvae are produced in yearly pulses and also included two additional life stages: the larval pre-competency and the non-reproducing juvenile stages (Gaylord and Gaines, 2000). The importance of better representing the species life stages is acknowledge recently in Bald et al. (2006) and Ramzi et al. (2006). The former proposes a system dynamic model for modelling the population dynamics of the gooseneck barnacle population in order to analyse different management strategies to the sustainable exploitation of resource in a marine reserve. The barnacle population is divided into three classes (juveniles, non-exploited adults and exploited adults) and eight different cases (management decisions) are simulated. The authors point out some model assumptions which could be subject to improvement. They mainly emphasize the high impact of the planktonic phase on the ulterior juvenile and adult structured and propose future improvement through the inclusion of environmental factors affecting recruitment and biomass. In Ramzi et al. (2006), a spatiotemporal mathematical model is developed to understand the larval ecology of the Moroccan sardine. It is an age-structured model in which two larval phases are considered (endogenous and exogenous) where processes such as transport (diffusive and advective), feeding and recruitment are involved. The transport parameters are assumed constant. Although no numerical simulation

27

is performed, the authors highlight the need of coupling the larval transport with a circulation model for recruitment estimates. They also suggest to extend the model to the exploitable phase to cover the whole sardine life cycle. In this work we propose a step towards a general stagestructured integrate model to understand the interplay between physical and biological processes of the population dynamics of organisms with complex life cycle in which the adult sessile phase is preceded by a planktonic larval phase. It extends Gaylord and Gaines (2000) although dealing with more details compared to (Alexander and Roughgarden, 1996; Gaylord and Gaines, 2000; Ramzi et al., 2006), taking into account the quality of the settlement habitat, a more flexible numerical model and a diagnostic model to recover the velocity field when the flow behavior is known at some discrete values in the domain of interest and at a given time. The numerical approach used in this work is based on the Finite Element Method (Hirsch, 1988; Hughes, 2000). The models (Alexander and Roughgarden, 1996; Gaylord and Gaines, 2000; Possingham and Roughgarden, 1990) were all numerically solved using the Finite Difference Method. The first step to solve approximately a given differential equation using the Finite Element Method is to reformulate the given problem (strong formulation) into a variational equivalent one (weak formulation) (Rektorys, 1980). Afterwards, based on the variational formulation, accuracy and stability requirements, the finite dimensional spaces are built using piecewise polynomials. Besides having a solid theoretical foundation, the finite element method can handle complicated geometries, general boundary conditions and variable or non-linear properties in a relatively easy way as compared with finite difference method (Hughes, 2000; Johnson, 2000). When dealing with real problems which typically have irregular geometry and both natural and essential boundary conditions, these properties are specially advantageous, avoiding additional geometric mapping and fictitious nodes as used in (Gaylord and Gaines, 2000). Thus, the proposed methodology allows flexibility in dealing with complex domains and is easily combined with an adaptive mesh procedure. The adaptive finite element method can be used to simulate the change in the domain when upwelling and relaxation events are present (Rio Doce et al., 2005) and to save computational efforts by adapting an unstructured finite element mesh to improve accuracy. The fully discrete problem is built based on the Galerkin method in space using piecewise bilinear finite element interpolation functions, and the implicit Euler method to approximate the time derivative. The discrete larvae transport equation, and the adult dynamics form a set of nonlinear algebraic equations. We propose here a predictor-corrector procedure to solve this problem. More sophisticated variational formulations can be easily introduced to deal with advection-dominated phenomenon ˜ et al., 2004). (Galeao The physical model is a diagnostic one (Montero and San´ın, 2001) that represents the velocity field based on the mass conservation principle. The velocity field is reconstructed throughout discrete velocity values in the domain of interest and at a given time. The proposed methodology is built following two steps. In the first one a simple interpolation method is used to determine an initial solution for the diagnostic model (see Kratzer et al., 2006) for a general study of

28

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

ment of the larvae (Fig. 1(b)). The larvae are produced by the reproducing adult population; they disperse in the water column by turbulent diffusion and advection (characterized by the velocity field v ≡ (vx , vy )) and they settle down when they strike the coast. The model assumes that the larvae remain homogeneously distributed in a water column of constant depth. The biological parameters for both larvae and adults are age-integrated averages (Possingham and Roughgarden, 1990). Hence, the density-independent larvae mortality rate is  and along the coast adults die (natural mortality and predation) at a rate  and produce larvae at a rate m. Larvae settle on the coast at a rate that is proportional to the product of the amount of free space available at time t, defined as F(y, t), and the larval concentration in the water column immediately adjacent to the shoreline at time t, L(0, y, t). The constant of proportionality that quantifies the amount of larvae to settle onto the shore is c. The space constraint is density dependent and can be represented by

Fig. 1 – Schematic representation of the model: (I) suitable habitat; (II) unsuitable habitat.

interpolating methods for current velocity). This velocity field is then variationally reconstructed based on the mass conservation principle using the diagnostic model. This methodology uses a steady data set. However, it can be straightforwardly applied to a time series through which flow changes would be more realistically represented. In such case, an evaluation of the time scales in which the physical and biological processes happen would be required. The outline of this paper is as follows. In Section 2 the mathematical model is presented, including the biological and the diagnostic models. Section 3 describes the proposed numerical methodology. Numerical experiments are conducted in Section 4 to show the robustness of the proposed numerical methodology. Finally, a discussion is drawn in Section 5.

2.

Mathematical model

2.1.

Population dynamics model

In this section, a mathematical model is briefly presented to describe the population dynamics of a single species in a twodimensional space (x ≡ (x, y) ∈ 2 ) and in time (t ∈ (0, T]). The organism has two main life phases: a sessile adult phase (B ≡ B(y, t)) restricted to the shore, and a larval phase (L ≡ L(x, y, t)) influenced by the pattern of ocean circulation (Alexander and Roughgarden, 1996). The main assumptions of the model are described as follow. The coast is described by a straight line and coincides with axis y. Initially, it is considered that only the central part of the coast is formed by an appropriate substratum for the settlement of larvae, as used in (Gaylord and Gaines, 2000; Possingham and Roughgarden, 1990), and that the population of adults is located along this habitat, as shown in Fig. 1(a). Later, the central part of the coast will be divided in two suitable habitats (I) and one not suitable (II) for the settle-

F(y, t) = A(y) − aB(y, t),

(1)

where A(y) is the total available free space for recruitment at position y and aB(y, t) quantifies the amount of occupied space at position y and time t. In this expression, a is the average basal area of an adult individual. The rate of change of the number of adult barnacles on the coast (position y) with respect to time is a balance between barnacle birth due to larvae settlement and adult mortality, as expressed by dB(y, t) = cF(y, t)L(0, y, t) − B(y, t). dt

(2)

It is assumed that the adult population on the coast is known at the initial time so as B(y, t = 0) ≡ B0 (y) represents the initial condition for the ordinary differential Eq. (2). The larvae produced by the adults are considered passive particles in the ocean and suffer dispersion through turbulent diffusion. The diffusion coefficient k (or turbulent diffusion coefficient) is assumed constant. The larvae can also be transported by advection due to currents in the water column; this movement is defined by the velocity field, the components of which are vx (velocity component in x direction) and vy (velocity component in y direction). With these assumptions, the rate of change of the larval population, L(x, y, t), with respect to time is described by the following equation: ∂L(x, y, t) ∂L(x, y, t) ∂L(x, y, t) + vx (x, y, t) + vy (x, y, t) = ∂t ∂x ∂y

 =k

∂2 L(x, y, t) ∂2 L(x, y, t) + 2 ∂x ∂y2

 − L(x, y, t) ,

(3)

x ∈ (0, xf ), y ∈ (0, yf ), t > 0, where xf is the frontal boundary that acts as a convergence zone where barnacle larvae accumulate (Alexander and Roughgarden, 1996). We assumed that the initial distribution of larvae is L(x, y, t = 0) = L0 (x, y).

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

29

Some hypotheses are assumed for the establishment of the boundary conditions. Along suitable habitats on the coast, larvae settle and are produced by the adult stage of the population. The larval flux on the coast (in an outwardly normal direction) is due to the larval production mB(y, t) and the settlement of larvae cL(0, y, t)F(y, t). Thus, the boundary condition that specifies the rate with which larvae leave and enter the water column on the coast is

 vx (x, y, t)L(x, y, t) − k

   

∂L(x, y, t) ∂x

x=0

= m B(y, t) − cF(y, t)L(x = 0, y, t).

(4)

It is assumed that the boundary at x = xf is a convergence zone where barnacle larvae accumulate so as the following null larval flux is imposed:

 vx (x, y, t)L(x, y, t) − k

∂L(x, y, t) ∂x

   

= 0.

(5)

x=xf

The diffusive larval flux at the south (y = 0) and north (y = yf ) boundaries are null and are expressed by



−k

∂L(x, y, t)   ∂y



=0 y=0

and

−k

∂L(x, y, t)   ∂y

Fig. 2 – Schematic representation of the stage-structured model (Gaylord and Gaines, 2000).

= 0.

(6)

y=yf

It should be remarked that the Eqs. (2) and (3) are coupled through the coastal boundary condition (4). The models (Alexander and Roughgarden, 1996; Possingham and Roughgarden, 1990; Roughgarden et al., 1988) assumed that larvae are produced continuously by adults and can settle immediately if they find appropriate substrate. Moreover, when they settle, they become adults immediately, and give birth to their own offspring instantaneously. These assumptions disregard some important features that affect the extent and pattern of dispersal of most marine species (Gaylord and Gaines, 2000) such as larval pre-competency and larval competency stages. The pre-competency phase is a period of larval development during which they are incapable of settling while the competency stage is a finite period of larval development during which larvae are able to settle into suitable habitat. According to (Gaylord and Gaines, 2000), the models (Alexander and Roughgarden, 1996; Possingham and Roughgarden, 1990; Roughgarden et al., 1988) underestimate the larval settlement and larval production rates and overestimate the adult mortality rate to compensate the absence of these life stages. Gaylord and Gaines (2000) argued that these modifications cannot replace a more realistic representation of life cycle of these organisms and developed a new model that includes two additional life stages. We compared these two models for the larval dispersion in a one-dimensional domain when upwelling and relaxation events are present. Although the two-phases model succeeded in qualitatively representing the adult population on the coast according to the flow changes, the four-phase model highlighted the importance of the reproduction lags and developmental periods on their distribution. The four-phase model proposed in Gaylord and Gaines (2000) is used in this work, as shown in Fig. 2. In the first

stage, larvae are produced in yearly pulses at the beginning of the reproductive season. This feature more accurately represents patterns observed in nature for many species, for which offspring are produced seasonally rather than continuously. The released larvae initially enter a pre-competency stage that lasts about 3 weeks. During this period larvae are transported by currents and dispersed by turbulent diffusion but cannot settle even if they encounter a suitable habitat. After that, larvae move to a competency stage, during which those larvae that reach suitable shoreline do settle, exiting the water column and joining the benthic population. This stage is finite and lasts as long as the pre-competency stage. Afterwards, all larvae that remain in the water column at the end of the competence period do die. The settled larvae pass through a several-month juvenile stage that lasts until the following reproductive season, when they release their offspring for the first time. The physical and biological parameters used in this model are described in Table 1(Gaylord and Gaines, 2000). The stage-structured model sets the larval production rate (adult fecundity) (m) and larval settlement coefficient (c) as functions of time. The larval production rate is a pulse at the beginning of each reproductive season and the larval settlement coefficient is zero except during the competency period. The eddy-diffusion coefficient (k) yields the dispersion of the larvae in the open ocean due to turbulence and depends on the phenomenon scale. Experiments following dye diffusing in the open ocean on the appropriate spatial and temporal scales for barnacle larvae estimate k = 10 cm2 /s. The reader may see Largier (2003) for a more precise scale analysis.

2.2.

Flow model

A diagnostic model is used for constructing a locally conservative velocity field. It is based on the mass conservation principle and assumes that velocity data are known from field measurements. The 3D diagnostic model for wind adjustment is presented, and Montero and San´ın (2001) is applied here for determining the 2D water velocity, assuming that there is no transient flow behavior. Under the constraint determined by the continuity equation, the adjustment is performed by

30

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

Table 1 – Physical and biological parameters Parameter

Nomenclature

Larvae per 100 m2 Adult per 100 m length of coastline Free space available for settlement per 100 m length of coastline Adult basal area Mortality rate of larvae Mortality rate of adults Larval production rate Larval settlement coefficient Larval pre-competency duration Larval competency duration

L(x, y, t) B(y, t) F(y, t) a   m(t) c(t) d1 d2

minimizing a least-square function. The Lagrange multiplier technique is used and the resolution of the saddle point problem yields an Euler–Lagrange equation that is solved using the finite element method. As the water column is assumed to have a constant depth, the initial field is obtained following a simple surface interpolation procedure as used. Assuming that the water density is constant and there is no time-dependent term, the continuity equation defined in the computational domain described in (3) is given by ∇v = 0

in ,

(7)

where v ≡ (vx , vy ) is the velocity field. The idea is to determine the velocity field that is closest to the field data, v0 ≡ (v0x , v0y ), and satisfies (7). The solution (vx , vy ) is obtained finding the saddle point (vx , vy , ϕ) of the functional given by E(vx , vy , ϕ) =

 

2

+ϕ

 ∂v

x

∂x

+

∂vy ∂y

 d,

(8)

the number of observation sites (stations). An interpolation technique is used to determine v0 (x, y) based on the square distances between the point (x, y) and the stations, denoted by rs , s = 1, . . . , OS. This procedure is called surface 1/r2 interpolation and follows the potential flow concept, being suitable for stable flows. Thus, after defining the domain discretization, the initial velocity field v0 = (v0x , v0y ) at position (x, y) is given by OS

v0 (x, y) =

and

vy = v0y +

1 ∂ϕ , 2 ∂y

3.



∂v0y ∂v0x + ∂x ∂y

(9)



(10)

To solve this problem it is necessary to append some boundary conditions. In general, imposing an homogeneous Dirichlet boundary condition (ϕ = 0) would lead to nonzero normal derivative, implying in mass flowing through boundary. Hence, defining n as the outward normal to the boundary, the following boundary conditions are used: ∇ϕ n = 0 ϕ=0

 +

 in .

(11)

. 1/rs2

Numerical method



where the Lagrange multiplier derivatives play the role of the velocity perturbation. Thus, the following elliptic equation in ϕ is obtained by substituting (9) into (7): ∂2 ϕ ∂2 ϕ + 2 = −2 2 ∂x ∂y

OS

s=1

+k 1 ∂ϕ 2 ∂x

vs /rs2

s =1

∂L  d − ∂t 

where ϕ is the Lagrange multiplier. This yields: vx = v0x +

1 × 10−4 m2 5.6 × 10−7 s−1 2.2 × 10−8 s−1 0 or 3.2 × 10−3 s−1 0 or 5 × 10−5 s−1 1.8 × 106 s 1.8 × 106 s

Defining S = H1 (), the usual Sobolev space of order 1, a variational formulation for (3) reads: for a given t ∈ (0, T] ≡ I and ∀ ∈ S, find L ∈ S such that

2

[(vx − v0x ) + (vy − v0y ) ]



Value and unit

on the coast (c ); otherwise (1 ),

where c ∪ 1 ≡ , the continuous boundary of . The initial velocity field is obtained from the discrete velocity data obtained through experimental measurements. Let these be denoted by vs = (vsx , vsy ), s = 1, . . . , OS, where OS is



∂ vx L d + ∂x 

∂L ∂ d + ∂x ∂x 



yf

L d − 

 vy 

∂L  d ∂y



∂L ∂ d ∂y ∂y 

(12)

[(mB − cFL)]|x=0 dy = 0. 0

This weak formulation is obtained by weighting the larval transport Eq. (3) and integrating over the whole domain. Afterwards, the second-order diffusive terms are integrated by parts and the x-component of the first-order advective term is integrated by parts in the x direction in order to introduce the natural boundary conditions (4)–(6), yielding (12). The discrete problem is now formulated based on the Galerkin method (Johnson, 2000) using the finite element method in space and the finite difference method in time. The finite element mesh is a partition of the domain  that consists of Nel finite elements e such that  =

Nel 

e and

e=1

Nel

e = ∅ (empty set). Let

e=1

Sh ⊂ S be the discrete counterpart of S defined as Sh = {h ∈ C0 (); h |e ∈ P1 (e )}, where P1 is the set of interpolation polynomials of degree less than or equal to 1 defined in each element e (Hughes, 2000).

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

The time discretization is built selecting 0 = t0 < t1 < . . . < T such that In ≡ (tn−1 , tn ) is a partition of I and tn = tn − tn−1 . For the sake of simplicity we set tn = t. Using the implicit Euler method to approximate the time derivative (Hirsch, h h h 1988), yields: find Lh n ∈ S , ∀ ∈ S and ∀In , such that



 h Lh n  d − t





− t



∂h d + t ∂x





h mB − cFLh n 

 vy



h ∂Lh n ∂ d + ∂x ∂x 

+k t

yf

vx Lh n



∂Lh n h  d ∂y



h ∂Lh n ∂ d ∂y ∂y 

 |x=0

0



h Lh n  d

+  t





Bn−1 (y) + t c A Ln (0, y) . [ t c a Ln (0, y) + t  + 1]

(16)

The solution of (16) depends on knowing the density of larvae on the coast, Ln (0, y). To solve the non-linear system of Eqs. (13) and (16) we developed the following predictor–corrector algorithm:

h Lh n−1  d

dy =

depends on knowing the distribution of the adult population on the coast (Bn ) at each time-step tn . The adult distribution Bn = Bn (y, t = tn ) is also obtained using the implicit Euler method in Eq. (2) for every position y of the mesh nodes which are lying on suitable coast habitat, yielding:

Bn (y) = (13)

31



Algorithm 1. Predictor–corrector.

In the nth time interval In , the finite element functions are defined as

Lh n (x, y, tn ) =

Nnp

Nj (x, y) Lh j,n ,

for x, y ∈ ;

(14)

j=1

h (x, y) =

Nnp

Nj (x, y) jh ,

for x, y ∈ ,

(15)

j=1

where Lh is the nodal value of Lh for the spatial node j at j,n time tn . jh is the value of the weighting function for node j and Nnp is the number of nodal points of the finite element mesh. Nj is the spatial interpolation function for node j and is assumed to be piecewise bilinear. Substituting the finite element functions (14) and (15) into (13) yields a system of algebraic equations having Nnp unknowns, the solution of which

This algorithm involves two phases to be performed in each time-step n, n = 1, . . . , nt, where nt = T/ t is the total number of time steps. The iterative procedure starts with an initial guess for the values Bn (y) and Ln (x, y) at every grid point (x, y), denoted by B0n and L0n , where the superscript stands for the number of the iterative step. These guessed values are plugged j into Eq. (16). Afterwards, the updated values Bn (y) are plugged j j into Eq. (13). Then, the values Bn (y) and Ln (x, y) are used as new guesses for the next iterative step (j + 1). The iterative procedure is repeated n c times, where n c stands for the number of iterations in the corrector loop. For transient problems, the convergence at each time is reached when the difference j j−1 between two subsequent solutions satisfies |Bn (y) − Bn (y)| < tol, where tol is the prescribed tolerance (Rio Doce et al., 2005). Finally, to determine the velocity field (9) to be used in (13) the elliptic problem (10) is solved for ϕ. To this end, the Galerkin finite element method is used. Defining the discrete

32

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

counterpart of the continuous space of test functions as V h = { h ∈ C0 ();

h |e ∈ P1 (e ),

h |1 = 0},

the discrete problem is: find ϕh ∈ V h such that

  

∂ϕh ∂ h ∂ϕh ∂ h + ∂x ∂x ∂y ∂y





 d = −

2 

∂v0y ∂v0x + ∂x ∂y



h d,

∀ h ∈ V h .

(17)

Using the finite element functions: ϕh (x, y) =

Nnp

Nj (x, y)ϕjh

and

h (x, y) =

j=1

Nnp

Nj (x, y) h j ,

j=1

in (17), a system of algebraic equations is obtained which has Nnp unknowns ϕjh . As we are not worried with computational efficiency, the algebraic equation systems are solved using a simple LU solver.

4.

Fig. 3 – Flow field I.

Numerical results

In this section, we present the numerical results obtained by applying the proposed methodology. In all experiments we consider that no larvae are in the ocean initially and 100,000 reproducing adults/100 m of coastline are homogeneously distributed along suitable coast habitats (10% cover of adults). Moreover, we use tol = 10−3 as the prescribed tolerance, t = 8 h as the time step-size, a uniform finite element mesh with 60 × 60 elements. We conducted a subset of runs with a variety of time step-sizes and finite element meshes in order to verify the convergence and stability of the proposed numerical procedure. This study is presented in Appendix A for the sake of clarity. To investigate the influence of ocean circulation and habitat quality on adult shoreline abundance and distribution of the barnacle Balanus glandula, we consider three different mass conservative flow field representations obtained from 12 discrete velocity field data. The positions of the stations and corresponding discrete velocity vector are given in Table 2. The role of the flow speed is also investigated by defining a characteristic velocity speed (v) for each data set. For each flow

Fig. 4 – Flow field II.

field, we conduct a series of runs using v equal to 0.3, 1 and 3 cm/s. The first flow field is shown in Fig. 3 and represents a flow with an eddy-circulation. The dot points indicate the posi-

Table 2 – Position of the stations (km) and their current velocities (cm/s) Station

Field I

Field II

s

x

y

vsx

1 2 3 4 5 6 7 8 9 10 11 12

1 6 1 8 1 5 9 1 15.4 9.6 1 9

95 88 72 70 64 52 42 28 28 16 7 7

0 0 0 0 0 v/2 v 0 0 −v −v/10 −v/10

vsy

x

−v −v −v −v −v −v 0 v −v 0 0 0

8 1 1 11 3 1 11 1 3 1 11 8

Field III

y

vsx

vsy

x

y

vsx

95 90 68 64 55 48 48 42 41 28 10 1

0 0 0 0 0 −v −v 0 0 0 0 0

v v v v v 0 0 −v −v −v −v −v

12 1 1 4 1 5 6 1 1 3 2 14

88 84 72 72 64 60 60 56 54 40 30 8

0 0 0 0 0 0 0 0 0 v v v

vsy −v −v −v −v −v −v −v −v −v 0 0 0

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

33

Fig. 5 – Flow field III.

tion of the stations. The second flow field situation is depicted in Fig. 4. It simulates a current flow that strikes the center of the coast in a perpendicular way (divergent flow). In the third flow field representation shown in Fig. 5, the flow has an L shape, which means that the flow closest to the coast has an alongshore southward movement until the middle of

Fig. 7 – Flow field I (0.3 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

the coast where it moves offshore. According to Gaylord and Gaines (2000), these three flow patterns simulate dominant features of several common large-scale flow fields. Each flow pattern is simulated considering the following coast conditions: a continuous habitat (Fig. 1(a)), where 50% of the coastline central part is formed by a suitable habitat; a sectioned habitat (Fig. 1(b)), such that the central coastline is divided into three sections, two suitable for settling (sections I - north and south) and one unsuitable (section II). In the following simulations, the adult population distributions refer to the end of the juvenile stage for each considered year.

4.1.

Fig. 6 – Flow field I (1 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

Flow field I

Fig. 6(a and b) illustrates the adult shoreline distribution resulting from the flow field I and v = 1 cm/s in the continuous suitable habitat and in the sectioned habitat, respectively. For the continuous habitat, Fig. 6a shows the growth of the adult population towards an equilibrium level, reached after about 10 years. The species concentrate primarily in the southern part of the suitable habitat, adjacent to the eddy, and the population decreases to extinction (defined as a percentage of adult cover below 0.02% Gaylord and Gaines, 2000) in the northern part. This occurs due to the larval transport. Ini-

34

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

and 1.4 × 10−3 , respectively, for the continuous habitat, and 1.6 × 10−1 , 7.9 × 10−3 and 4.0 × 10−4 , for the sectioned habitat. These results strongly suggest that the equilibrium state eventually will be reached.

4.2.

Flow field II

Fig. 9(a and b) presents the results with the flow field II and v = 1 cm/s. For the continuous habitat, Fig. 9(a) shows that the adult population increases very fast from the initial condition (10% cover) to an equilibrium level (approximately 55% cover), reached in 5 years. For the sectioned habitat, we can observe a very different result (Fig. 9(b)). In both southern and northern sections, the adult population slides southward and northward, respectively, before going to extinction. This indicates that the persistence of the barnacle population is due to the adults present in the central region, adjacent to the perpendicular current, where larvae are transported against the coast. In the first case, the larvae produced by the central adult population during the reproductive season can be swept southward and northward when the current splits. In this way, they can settle when in contact with the remaining shore. In the second case, all the larvae produced in sections I are transported away from suitable habitat.

Fig. 8 – Flow field I (3 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

tially, the larvae are transported downstream and, later, they fall within the recirculation zone. As the flow is assumed steady, they can return to the shore after the end of the precompetency period. For the sectioned habitat, Fig. 6b shows a similar result. In the southern section, the adult population increases until an equilibrium level is reached while the adult population goes towards extinction in the northern section. As the velocity speed decreases, the adult shoreline distribution tends to be velocity-independent, as shown in Fig. 7(a and b) for v = 0.3 cm/s. In these cases, the adult abundance is almost uniform and this is primarily due to the habitat quality. On the other hand, when v increases to 3.0 cm/s, the species concentrate primarily in the southern coast region, adjacent to the eddy (Fig. 8(a and b)), as similarly happened for v¯ = 1.0 cm/s. In this case, the simulation time reached year 90 and, after an initial growth, the southern region adult population smoothly declines towards an equilibrium state. Although it was not get at year 90, the maximum relative difference (among the coastal nodal points) between adult population density at these two subsequent years 29–30, 59–60 and 89–90 are 4.2 × 10−1 , 1.9 × 10−2

Fig. 9 – Flow field II (1 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

35

A decrease in the velocity speed has similar consequences for the flow pattern III, as shown in Fig. 13(a and b) for v = 0.3 cm/s. Because of the flow pattern, larvae are swept faster from the coast when v = 3.0 cm/s so as the adult population declines rapidly to extinction for both the continuous and the sectioned habitats (Fig. 14(a and b)). The previous results agree qualitatively well with the ones presented in Gaylord and Gaines (2000), and highlight the great impact of circulation pattern on population dynamics of a species that has a dispersing larval stage. Depending on each flow pattern, there exist critical current speeds for population persistence. Thus, the persistence and abundance of the benthic population do not depend exclusively on the habitat quality.

5.

Discussion

In this work we present a simple and general finite element methodology to solve the population dynamics of a marine organism (barnacle B. glandula) with a complex life cycle. Our previous works (Rio Doce et al., 2005) together with the convergence tests presented in Appendix A indicate that the numerical results are accurate and stable.

Fig. 10 – Flow field II (0.3 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

As for flow field I, a decrease in the velocity speed also decreases the dependence of the adult shoreline distribution on the flow pattern, as shown in Fig. 10(a and b) for v = 0.3 cm/s. The adult abundance is more uniform than in the previous flow pattern and depends almost exclusively on the habitat quality. When velocity speed increases to 3.0 cm/s the species either reaches an equilibrium state for the continuous habitat (Fig. 11(a)) or goes toward extinction for the sectioned habitat (Fig. 11(b)) in a very fast way.

4.3.

Flow field III

Fig. 12(a and b) shows the adult shoreline distribution resulting from the flow field III. The released larvae are swept southward until the middle of the coast and, after that, offshore. For both the continuous suitable habitat and the sectioned habitat, the adult population declines to extinction. In the first case (Fig. 12(a)), the biggest adult abundance occurs in the central region of the habitat at the very beginning. In the other case (Fig. 12(b)), the population slides downstream in the northern section while in the southern section the population declines to extinction more uniformly.

Fig. 11 – Flow field II (3.0 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

36

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

Fig. 12 – Flow field III (1 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

This species has a planktonic larval stage that is dispersed due to turbulent diffusion and advection. These physical processes strongly influence the adult abundance and distribution on the coast (Roughgarden et al., 1988; Possingham and Roughgarden, 1990; Alexander and Roughgarden, 1996; Gaylord and Gaines, 2000) and thereby affecting other species dynamics with which the adult population interacts in the benthic habitat (Connolly and Roughgarden, 1998, 1999). Due to the critical role that oceanographic processes play on the populations dynamics and community structure, the present study incorporates a simple model to better represent coastal currents. Therefore, its effects on the larval dispersal patterns of rocky intertidal organisms with planktonic stages may also be better represented. The used approach allows the construction of physically relevant current patterns which otherwise could not be reproduced without using an oceanographic model Ramzi et al. (2006). In principle, any flow pattern defined by discrete velocity values in the domain can be physically and consistently represented and may be easily coupled to the biological model. We present several scenarios that allow the understanding about the qualitative effects of idealized flow patterns and advection levels over the population dynamics on the coast. They essentially show that coastal flows may affect

Fig. 13 – Flow field III (0.3 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

the adult distribution independently of the habitat quality. This work is a step towards applying the developed methodology to more complex scenarios that can assimilate real unsteady velocity fields in more complex geometries (Graham and Largier, 1997; Largier, 2003). The developed approach allows the introduction of many other modeling changes in a quite straightforward way so as to qualitatively evaluate the interplay between biological and physical processes. The model can be evaluated under a vast range of parameters, non-uniform environment and parameter dependence on some other physical property, such as temperature, for example. Also, different space and time scales may be considered (Gaines et al., 2003). Such qualitative analysis provides understanding about the population dynamics model and may yield valuable insights for constructing more general/realistic models. One actual application is the definition of the spatial design and/or study of marine-protected areas (MPAs), that are recognized among specialists to be important elements of marine conservation (Largier, 2003). Another important issue in marine community studies in a biogeographical context is the role of benthic processes such as competition and predation. Models of interacting

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

37

Appendix A

Fig. 14 – Flow field III (3.0 cm/s): adult shoreline distribution in the continuous habitat (a) and in the sectioned habitat (b).

species may be easily incorporated in the proposed numerical approach just yielding an increase of computational cost due to a larger non-linear system of equations. To overcome this drawback, a more efficient algorithm/solver may be introduced and there is a vast literature on this subject (Johnson, 2000). Thus, the main value of the proposed numerical methodology is to provide a reliable tool to study the distribution and abundance of sessile organisms with a planktonic larval stage.

In this section, we verify the convergence and stability of the proposed numerical solutions by conducting a subset of runs with a variety of time step-sizes and finite element meshes. Two numerical examples are considered. In the first one, different time step-size values ( t = 8, 4, 2, 1.11, 0.55 and 0.028 h) are used to solve a population dynamics of the barnacle B. glandula in a simple one-dimensional domain  = [0, xf ], discretized using 200 finite elements. We consider that the larval transport occurs by diffusion and advection (vx = 0.3 cm/s) and the front (xf ) is located at 19.2 km offshore. Fig. A.1 shows the variation of the percent of cover of reproducing adults on the coast with time. It shows that the solutions for decreasing time steps are convergent in time. They are slightly different only during the first 2 years and, after that, they reach a steady behavior characterized by the limit cycle. This cycle is characterized by an increase of the adult population due to settlement and a decrease of adult population due to adult and larval mortality. In the second example, the convergence of the finite element solution is studied using three different uniform meshes (4 × 30, 8 × 60 and 16 × 120 elements) in the two-dimensional domain presented in Fig. 1. In this case, we also consider the transport of larvae by diffusion and advection (alongshore flow at a constant velocity vy = 0, 3 cm/s) and t = 8 h. The results are depicted in Fig. A.2(a) to compare the effect of the mesh size on the solutions. The convergence is clearly reached when the mesh is refined. The same behavior is obtained for non-uniform meshes. This is illustrated in the solutions depicted in Fig. A.2(b) obtained keeping constant the number of finite elements in the y-direction and changing the number of elements in the x-direction (15 × 60, 30 × 60, 60 × 60 and 90 × 60 elements). We can observe that the finite element solutions on the coast for this flow pattern does not depend on the number of elements in the x-direction. For all the flow patterns and velocity speeds considered in this work, the numerical solution were stable, free of spurious modes. The reader may report to Rio Doce et al. (2005) for accuracy results.

Acknowledgments The first author would like to thank the fellowship pro´ vided by FAPERJ and project GEOMA—“Rede Tematica ˆ de Pesquisa em Modelagem da Amazonia”. The second author thanks the Brazilian Government, through the Agency CNPq, for the financial support provided. This work was also supported by the project PRONEX/FAPERJ E26/171.199/2003.

Fig. A.1 – Temporal convergence test.

38

e c o l o g i c a l m o d e l l i n g 2 1 4 ( 2 0 0 8 ) 26–38

Fig. A.2 – Convergence test for mesh refinement for variation of element in x-direction and y-direction (a) and in x-direction (b).

references

Alexander, S., Roughgarden, J., 1996. Larval transport and population dynamics of intertidal barnacle: a coupled benthic/oceanic model. Ecological Monographs 66 (3), 259–275. Bald, J., Borja, A., Muxika, I., 2006. A system dynamics model for the management of the gooseneck barnacle (pollicipes pollicipes) in the marine reserve of Gaztelugatxe (Northern Spain). Ecological Modelling 194, 306–315. Connolly, S., Roughgarden, J., 1998. Theory of marine communities: competition, predation, and recruitment-dependent interaction strength. The American Naturalist 151 (4), 311–326.

Connolly, S., Roughgarden, J., 1999. Theory of marine communities: competition, predation, and recruitment-dependent interaction strength. Ecological Monographs 69 (3), 277–296. Gaines, S., Gaylord, B., Largier, J.L., 2003. Avoiding current oversights in marine reserve design. Ecological Applications 13, S32–S46. ˜ A.C., Almeida, R.C., Malta, S.M.C., Loula, A.F.D., 2004. Galeao, Finite element analysis of convection dominated reaction diffusion problems. Applied Numerical Mathematics 48, 205–222. Gaylord, B., Gaines, S., 2000. Temperature or transport? Range limits in marine species mediated solely by flow. The American Naturalist 155, 769–789. Graham, W.M., Largier, J.L., 1997. Upwelling shadows as nearshore retention sites: the example of northern Monterey Bay. Continental Shelf Research 17, 509–532. Hirsch, C., 1988. Numerical Computation of Internal and External Flows. Vol. 1. John Wiley & Sons. Hughes, T.J.R., 2000. The Finite Element Method - Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc. Johnson, C., 1990. Numerical solution of partial differential equations by the finite element method. Cambridge University Press. Kratzer, J.F., Hayes, D.B., Thompson, B.E., 2006. Methods for interpolating stream width, depth and current velocity. Ecological Modelling 196, 256–264. Largier, J.L., 2003. Considerations in estimating larval dispersal distances from oceanographic data. Ecological Applications Supplement 13, S71–S89. Montero, G., San´ın, N., 2001. 3-D modelling of wind field adjustment using finite differences in a terrain conformal coordinate system. Journal of Wind Engineering and Industrial Aerodynamics 89, 471–488. Pineda, J., 2000. Linking larval settlement to larval transport: assumptions, potentials and pitfalls. Oceanography of the Eastern Pacific 1, 84–105. Possingham, H., Roughgarden, J., 1990. Spatial population dynamics of a marine organism with a complex life cycle. Ecology 71, 973–985. Ramzi, A., Hbid, M., Ettahiri, O., 2006. Larval dynamics and recruitment modelling of the Moroccan Atlantic coast sardine (sardina pilchardus). Ecological Modelling 197, 296–302. Rektorys, K., 1980. Variational Methods in Mathematics, Science and Engineering, 2nd ed. D. Reidel Publishing Company. Rio Doce, A.P.C., Almeida, R.C., da S. Costa, M.I., 2005. An adaptive finite element model to investigate the interplay between circulation patterns and the dynamics of a marine organism. Fourth Brazilian Symposium on Mathematical and Computational Biology/First International Symposium on ´ Mathematical and Computational Biology, Ilheus, Bahia. Vol. I. pp. 56–76. Roughgarden, J., Gaines, S., Possingham, H., 1988. Recruitment dynamics in complex life cycles. Science 241, 1460–1466. Roughgarden, J., Iwasa, Y., Baxter, C., 1985. Demographic theory for an open marine population with space-limited recruitment. Ecology 66 (1), 54–67.