Digital Commons @ Michigan Tech Digital Commons @ Michigan Tech A hybrid lagrangian–eulerian particle model for ecosystem A hybrid lagrangian–eulerian particle model for ecosystem simulation simulation

: Current numerical methods for simulating biophysical processes in aquatic environments are typically constructed in a grid-based Eulerian framework or as an individual-based model in a particle-based Lagrangian framework. Often, the biogeochemical processes and physical (hydrodynamic) processes occur at different time and space scales, and changes in biological processes do not affect the hydrodynamic conditions. Therefore, it is possible to develop an alternative strategy to grid-based approaches for linking hydrodynamic and biogeochemical models that can signiﬁcantly improve computational efﬁciency for this type of linked biophysical model. In this work, we utilize a new technique that links hydrodynamic effects and biological processes through a property-carrying particle model (PCPM) in a Lagrangian/Eulerian framework. The model is tested in idealized cases and its utility is demonstrated in a practical application to Sandusky Bay. Results show the integration of Lagrangian and Eulerian approaches allows for a natural coupling of mass transport (represented by particle movements and random walk) and biological processes in water columns which is described by a nutrient-phytoplankton-zooplankton-detritus (NPZD) biological model. This method is far more efﬁcient than traditional tracer-based Eulerian biophysical models for 3-D simulation, particularly for a large domain and/or ensemble simulations.


Introduction
Current numerical methods for simulating biogeochemical processes in aquatic environments are typically constructed as individual-based models (IBMs) using Lagrangian particles [1][2][3] or as grid-based concentration models in an Eulerian framework [4][5][6]. IBMs allow for a mechanistic description of individuals and of interactions among individuals, represented by an ensemble of particles. Each individual contains a set of state variables (e.g., age, size, and nutrient quota) with corresponding physiological traits and behavioral traits [3], and the population-level properties emerge as a result of the cumulative behavior of the individuals [7]. IBMs have rapidly gained popularity in ecological modeling, particularly when simulating complete life cycles, adaptive behavior, and intrapopulation variability in response to internal and external environmental conditions becomes essential [8,9]. Large computational demand has long been known as one of the hallmark problems for IBM [2,7,8]. This is still true even with increased computational resources and implementation of the strategic approaches for reducing the number of individuals explicitly simulated in the model such as introducing the super individual particles or representative spaces [2].
On the other hand, the traditional grid-based Eulerian approach is also widely used in coupled physical-biogeological modeling [6,7,10,11] in the inland water and ocean modeling communities [12][13][14][15][16][17][18]. Unlike the particle models in which mass is transported as discrete particles [19,20], the Eulerian (concentration) model assumes average properties (state variables) of a population within a control volume, and estimates the change of the property using mass balance and reaction equations [21]. It does not describe intrapopulation variability but focuses on characteristics of the population mean, which is appropriate when growth kinetics are formulated as a function of the external conditions [7].
In general, the Eulerian model employs finite difference or finite element schemes to solve the governing reaction-transport equations [4,6]. Equations for the time evolution of state variables of the biophysical model include advection and diffusion terms which depend on hydrodynamic variables, as well as source and sink terms representing growth, decay, and interaction with other biogeochemical variables [6,11,22]. The property concentration fields (C i , i = 1, 2, 3, . . . .) are often calculated using a set of advection-diffusion equations: where D is the total water depth, u, v, and w are the x, y, and z components of the water velocity, K h is the vertical thermal diffusion coefficient, F c is the horizontal diffusion term, and C i,source and C i,sink represents the sources and sinks of C i , respectively, due to the biological processes which are typically described using a set of biological process equations. However, even in the relatively simple, less computationally demanding Eulerian model, a major practical challenge is that the biological submodel often involves a large group of parameters for calibration and confirmation which requires a considerable amount of computational time to tune the model. As shown in Equation (1), tuning the simulation of biological processes (e.g., changes in parameterization, initial and boundary conditions) requires a complete time integration of the entire equation, although bio-physical coupled models can have different time steps for physical and biological process terms as they may have different time and space scales. However, biophysical processes are generally not two-way coupled. In other words, one can often assume that changes in biological processes (in our case, the resulting changes in nutrient-phytoplankton-zooplankton-detritus (NPZD) property concentration) do not affect the hydrodynamic condition (currents, temperature, mixing, etc.). This indicates that there may be a more computationally efficient approach to resolve the impact of hydrodynamics on the biological processes rather than explicit integration of Equation (1) for each biological component every time the biological submodel is tuned. The property-carrying particle model (PCPM) is developed to test the feasibility of an alternative strategy to grid-based approaches for linking hydrodynamic and biogeochemical models in the Eulerian framework that may reduce the problems mentioned above, particularly in regions where advection-diffusion plays a key role in regulating biogeochemcal processes. Instead of grid-based, time-averaging of hydrodynamic variables, the hydrodynamic model is used to calculate the Lagrangian trajectories of a large number of current-following tracer particles; these trajectories become the linking mechanism between the hydrodynamic model and the biogeochemical model. In the hybrid Lagrangian-Eulerian PCPM, each current-following tracer particle carries with it a number of time-varying properties which correspond to the state variables of the Eulerian biogeochemical model. The PCPM also employs its own horizontal grid system or series of regions which is independent of the hydrodynamic model grid and is used to calculate local average values of the particle-based properties. These cell-based properties allow all particles within a PCPM cell to influence the properties of other particles within the same cell or region and allow for display and analysis of biogeochemical fields. PCPM also differs from typical IBMs in that the tracer particles in PCPM typically carry 'field' based properties like concentration, as opposed to properties associated with an individual organism.
PCPM uses a computational grid system which is independent of the grid system used to compute currents for particle trajectories. The PCPM computational cells are used to define regions in which the properties carried by the particles are allowed to interact with one another. In this respect, PCPM is conceptually similar to the classic particle-in-cell (PIC) method with PM (particle-mesh) interactions [23,24], popularly used in plasma simulation [24,25]. PIC methods can also be mesh-independent by allowing direct particle-particle (PP) interactions, or a combination of PM and PP [23][24][25][26]. In PCPM, a basic simplifying assumption is that only particles within a single PCPM cell are allowed to interact, such as the PIC-PM method. The advantage of this approach is that it is conceptually intuitive to implement and computationally efficient to program.
An alternative approach to implementing a PCPM would allow particle-based properties to influence particle trajectories, perhaps through buoyancy or sinking. In this case, the PCPM would have to be coupled directly with the particle trajectory calculation. In the initial implementation of PCPM in this paper, we consider only the uncoupled case.
To illustrate more clearly the type of application envisioned for PCPM, we constructed and applied a rudimentary biophysical model to an actual aquatic system, Sandusky Bay, where the physical transport of flow and nutrient loading from the Sandusky River has proven to be critical to the ecosystem [27,28]. Since the mid-1990s, harmful algal blooms (HABs) have become the new norm for summer months in the Lake Erie ecosystem [29]. Harmful algal blooms occur in the system when cyanobacteria are provided the right temperature, light, and nutrient conditions to proliferate [30]. When these blooms transpire, they have many adverse impacts. At the local ecosystem level, HABs result in depleted dissolved oxygen levels below the lake's surface threatening the survival of organisms living below the surface [31]. Additionally, some cyanobacteria species produce a toxin, such as microcystin, which affects the nervous system, liver, and kidney further impeding aquatic organisms and humans [31].
The colonial cyanobacterium Microcystis dominates the cyanobacterial community in the offshore waters of western Lake Erie; however, the filamentous cyanobacterium Planktothrix dominates the nearshore bays and tributaries [28], such as in Sandusky Bay [27,28,32]. Situated on Lake Erie's south-western coast, Sandusky Bay borders Ohio's Ottawa, Erie, and Sandusky counties (Figure 1). From a physical aspect, Sandusky Bay is relatively shallow bay with an average depth of roughly 2.6 m as well as occupying a relatively small area [32]. The primary draining watershed to Sandusky Bay originates from the Sandusky River at the west end of the bay. The Sandusky River drains a 1420 square mile area; of which, over 80% is dedicated to agricultural production [31]. This largely agricultural watershed leads to high nitrogen and phosphorus entering Sandusky Bay. Combining these high nutrient loads with the physical aspects leads to very high concentrations of nitrogen and phosphorus within Sandusky Bay, thus resulting in these cyanobacteria blooms (Planktothrix agardhii) [27,28,32].
The intent of this study is to test the feasibility of PCPM for biological-physical coupled modeling and examine its effectiveness and computational efficiency in practical application by implementing it in relation to HABs in Sandusky Bay. The remaining sections of this paper are organized as follows: details of PCPM and the design of numerical experiments are described in Section 2. The model results of idealized cases and the practical application to Sandusky Bay are presented in Section 3. A discussion and summary of the PCPM applications are included in Section 4. Figure 1. Sandusky Bay is situated on Lake Erie's south-west coast occupying a small portion of the Great Lake's coastline. Sandusky Bay is relatively shallow bay with an average depth of ~2.6 m. The primary draining watershed to Sandusky Bay originates from the Sandusky River on the west end of the bay. Sampling stations ODNR1 and EC1163 are denoted with green dots.

Property-Carrying Particle Model (PCPM)
In this implementation of PCPM, particle trajectories are pre-computed based on the output of a hydrodynamic model and are independent of the particle properties. An initial particle density (i.e., total number of particles/volume of computational domain) is selected and particles are randomly distributed throughout the computational domain. Particles are not allowed to leave the computational domain except at hydrodynamic outflows. At hydrodynamic inflows, new particles are introduced with the same density as the initial distribution. The total number of active particles is not strictly preserved, but if there is a net balance of hydrodynamic inflows and outflows, the total number of particles is approximately constant.
Any suitable method can be used to generate the Lagrangian particle trajectories. Typically, the trajectories are calculated from a time integration of the Lagrangian equations of motion: where (x, y, z) is the particle's position in 3 dimensions, (u, v, w) is the local fluid velocity vector, and t is time. For the two idealized examples presented in this paper, the trajectories are calculated semianalytically from a simple, idealized flow field. The third, more realistic example, demonstrates the use of a full hydrodynamic model of a natural basin (i.e., Sandusky Bay) to compute currents and trajectories.

Property-Carrying Particle Model (PCPM)
In this implementation of PCPM, particle trajectories are pre-computed based on the output of a hydrodynamic model and are independent of the particle properties. An initial particle density (i.e., total number of particles/volume of computational domain) is selected and particles are randomly distributed throughout the computational domain. Particles are not allowed to leave the computational domain except at hydrodynamic outflows. At hydrodynamic inflows, new particles are introduced with the same density as the initial distribution. The total number of active particles is not strictly preserved, but if there is a net balance of hydrodynamic inflows and outflows, the total number of particles is approximately constant.
Any suitable method can be used to generate the Lagrangian particle trajectories. Typically, the trajectories are calculated from a time integration of the Lagrangian equations of motion: where (x, y, z) is the particle's position in 3 dimensions, (u, v, w) is the local fluid velocity vector, and t is time. For the two idealized examples presented in this paper, the trajectories are calculated 1.
Read particle locations (x, y, z) and temperature. This step simply updates the location of each particle that is being used in the computation. Figure 2 is a conceptual representation of a PCPM computational cell, Particles (m 1 , m 2 , m 3 , . . . ) move in and out of the cell at each PCPM time step based on their trajectories as computed from the hydrodynamic model. The total number of particles for a particular computation is assumed to be fixed for the duration of the computation, although some particles may enter or leave the PCPM domain during the computation. Water temperature or other physical properties from the hydrodynamic calculation can be stored along with the pre-computed particle trajectories and can be included as one of the properties (P1, P2, P3, . . . ) carried by the particle.

2.
Determine the PCPM cell for each particle. In Figure 2, the PCPM cell is represented by the enclosing rectangle. The PCPM domain need not coincide with the domain that was used for the hydrodynamic simulation and computation of particle trajectories. It can be regular or irregular, as long as there is a prescribed method to calculate which PCPM cell contains a prescribed particle position (x, y, z). The PCPM cells are the volumes within which particle properties can interact, that is, during a single time step, all particles within a PCPM cell can influence the evolution of particle properties within that cell, but are independent of other cells.

3.
Apply boundary conditions to any particle-based properties that require them. If there is a property (e.g., concentration of a dissolved nutrient) that needs to be specified as a boundary condition, then particles within the cell where the boundary condition needs to be applied will have that property adjusted to meet the boundary condition. For example, in a cell that is associated with an inflow to the domain, the properties that are being carried into the domain through the inflow are adjusted to take account of the change in that property for particles within that cell. Alternatively, if particles from the hydrodynamic-based trajectory calculation are entering a PCPM cell, the values of the associated properties for each particle need to be specified.

4.
Calculate PCPM cell-based averages of each property. In this step, the averages of K th property for cell n are calculated as: where the summation includes all L particles (m 1 , m 2 ,...m L ) currently within cell n. L is the number of particles within that cell. If no particles are present in a particular cell, PCPM uses the values of PK n from the previous time step.

5.
Calculate the time evolution of the cell-based properties (and particle-based properties if necessary) using the process equation defined for that property. The process equations can incorporate terms which depend on either particle-based or cell-based properties, or both, i.e., Note that M indicates m 1 , m 2 ,...m L . The form of FN is completely general and depends on the problem being solved. For instance, in a NPZD model, the Pi, (i = 1, 2, 3, . . .) would be N, P, Z, D, and water temperature, and the FN would be the process equations relating these properties.
Since the cell-based averages have already been computed, the right-hand side of Equation (4) is independent of the left-hand side, so the computation of the evolution equations can be carried out in parallel. This is another key design feature of PCPM allowing it to take full advantage of multiprocessing computer environments, both symmetric multi-processing (SMP) and massively parallel processing (MPP).

6.
Redistribute cell-based properties to particles within each cell by replacing the particle-based property with a weighted average of the cell-based property. After the evolution equations have been carried out (Step 5), particles within an individual cell most likely carry a range of different values of the various properties, which vary around the new cell-based average computed in Step 5, PK n (t + ∆t). PCPM provides an optional step to reduce the variance of the new particle-based properties within each cell. This optional step is applied as a 'nudging' term, i.e., where 0 < α i < 1 is the redistribution weight (i.e., nudging) factor. If α i = 0, no adjustment is carried out and the particle-based property remains unchanged. If α i = 1, then all particles within a cell are assigned the cell-based average of that property. This step can be useful to smooth results when limited particle density results in excessive within-cell variability. Step 5, ( + ∆ ) . PCPM provides an optional step to reduce the variance of the new particlebased properties within each cell. This optional step is applied as a 'nudging' term, i.e., where 0 < < 1 is the redistribution weight (i.e., nudging) factor. If = 0, no adjustment is carried out and the particle-based property remains unchanged. If = 1, then all particles within a cell are assigned the cell-based average of that property. This step can be useful to smooth results when limited particle density results in excessive within-cell variability. Note that all steps except 3 and 5 are independent of the specific problem, i.e., they will be carried out the same way no matter how many properties are attached to the particles or what those properties represent. More importantly, steps 1 and 2 only need to be run once regardless of modifications in biological processes at the later stage. These are two of the key designs of PCPM that contribute to enhanced computational efficiency.

Idealized Case 1: Advection-Diffusion Plume
In PCPM, diffusion is provided mainly by particle trajectories, although the cell-based averaging of particle properties and the (optional) redistribution of cell-based properties to particles within the cell can also act as diffusive terms. To demonstrate the effect of particle trajectory diffusion on particle properties, we constructed a 500 m wide × 2000 m long channel divided into 10 m square cells. Note that all steps except 3 and 5 are independent of the specific problem, i.e., they will be carried out the same way no matter how many properties are attached to the particles or what those properties represent. More importantly, steps 1 and 2 only need to be run once regardless of modifications in biological processes at the later stage. These are two of the key designs of PCPM that contribute to enhanced computational efficiency.

Idealized Case 1: Advection-Diffusion Plume
In PCPM, diffusion is provided mainly by particle trajectories, although the cell-based averaging of particle properties and the (optional) redistribution of cell-based properties to particles within the cell can also act as diffusive terms. To demonstrate the effect of particle trajectory diffusion on particle properties, we constructed a 500 m wide × 2000 m long channel divided into 10 m square cells. Particles were introduced at random locations along the center 400 m section of the left edge of the channel at the rate of 100/s. The particles were assigned an along-channel velocity of 2 m/s. Horizontal diffusion was added using a random-walk perturbation to the particle trajectories of 2r √ 2k h ∆t in both cross-channel and long-channel directions. Here, r is a uniformly distributed random number in the range [−1, 1], k h is the horizontal diffusion coefficient (10 m 2 /s in this experiment), and ∆t is the time step for the particle trajectory calculation (1 s).
In this example, PCPM particles carry only one property, concentration (P1 = C), and there is no time evolution equation (step 5, above). The purpose of this example is to illustrate how PCPM simulates horizontal diffusion through a combination of the particle trajectories and the cell-based averaging in step 6. To simulate a concentration plume, particles introduced in the center of the left wall (−50 m < y < 50 m) are assigned the initial condition C = 1. Particles entering the channel outside this region have an initial condition of C = 0.

Idealized Case 2: Vertical Settling
Since this implementation of PCPM does not allow the properties carried by the particles to influence particle trajectories, the question arises of how to simulate the vertical transport of a property when the vertical transport depends on the property itself, such as sediment settling or biologically generated buoyancy. In PCPM, the answer is simply to solve the vertical transport at the PCPM cell-based Eulerian framework in step 5 as a traditional cell-based method. Interaction of particle properties with adjacent cell averages is technically not allowed in the basic PCPM framework, but an exception is made in this case. The vertical advection-diffusion equation for sediment concentration is shown below: where w s is the bulk settling velocity of the suspended material and k z is the vertical diffusion coefficient.
Since vertical diffusion is already included in the particle trajectories, PCPM only needs to consider the first term on the right-hand side of (7) to account for the additional vertical transport that depends on the property itself. To implement this term in the PCPM, the process equation for a particle carrying a property C m in vertical cell k looks like: where C k (t) is the average concentration in vertical cell k, C k−1 (t) is the average concentration in the next higher vertical cell, and ∆z is the spacing between the centers of the cells. For particles in the top cell (k = 0), we set: and for particles in the bottom cell (k = kmax), we set: As a test case, we examine the vertical settling in a one-dimensional water column of depth d with particles moving vertically only through vertical diffusion. Particles are initially distributed randomly in the column and then move with a random walk velocity of 2r √ 2k z ∆t where r is a uniformly distributed random number in the range [−1, 1] and k z is the vertical diffusion coefficient. Particles are not allowed to cross the surface or bottom boundaries. Thus, in this experiment, the number of particles is constant and the particles are always approximately uniformly distributed in the vertical due to vertical mixing.
For the experiment, we set C = 1 as the bottom boundary condition by assigning this value at the beginning of each time step to all particles in the lower half of the bottom cell. The initial condition in other cells is C = 0. For the test case, we set the number of particles to 1000, d = 20 m, k z = 10 −4 m 2 s −1 , and the redistribution parameter α = 0.1. Three runs were made with 5, 10, and 20 vertical cells respectively. The PCPM is integrated in time with ∆t = 1 h.

Sandusky Bay Model
The hydrodynamic model used in this study is FVCOM (finite volume community ocean model) [33]. FVCOM is an unstructured-grid, finite-volume, three-dimensional (3-D) primitive equation ocean model with a generalized, terrain-following coordinate system in the vertical and a triangular mesh in the horizontal. The unstructured grid can be designed to provide a customized variable resolution to both coastline and bathymetry. With the merits of ideal geometric fitting and local refinement of mesh resolution, FVCOM has been used in numerous applications to estuaries, coastal oceans, and the Great Lakes [34][35][36][37][38][39][40]. These characteristics make the model well suited for the study of Sandusky Bay.
Although this study focuses on Sandusky Bay, FVCOM is configured to simulate physical dynamics for all of Lake Erie including a high-resolution Sandusky Bay-FVCOM developed in this study, thus providing reliable representation of large-scale background circulation and the role of remote forcing impacting the water movement in the bay through the opening; additionally, this configuration avoids the impact of setting an artificial numerical boundary condition for our target region. The hydrodynamic model is well calibrated for the Lake Erie, based on the next-generation National Oceanic and Atmospheric Administration (NOAA) Lake Erie Operational Forecast System (LEOFS; see Kelley et al., [40] for detailed model validation), a real-time nowcast and forecast model that is built on the FVCOM. In the upgraded NOAA operational model for Lake Erie [40], the FVCOM model is developed with horizontal resolution ranging from 100 to 2500 m, and 21 uniform vertical sigma (terrain-following) layers for Lake Erie. The advantage of our model setting is that model resolution varies from 100-2500 m (coarse) in the open lake to 10-50 m (fine) in Sandusky Bay, affording a high degree of resolution across the 20 km × 3 km study site and adequately resolving the geographic complexity and coastal hydrodynamic conditions of that system (Figure 3). The model configuration yields a total of 73,000 grid elements (cells) in the horizontal plane with 50,000 of them resolving the bay.
In the PCPM implementation, 86,000 initial particles are randomly distributed throughout Sandusky Bay with a total water volume of 3.01 × 10 8 m 3 . With PCPM cell resolution of 200 m × 200 m and the mean water depth of 2.6 m in Sandusky Bay, each PCPM cell contains 30 particles on average. New particles are introduced from the Sandusky River with the same density as the initial distribution. The number of new particles released from the river mouth varies greatly in accordance with the river flow rate. Table 1 presents the number of new particles released each month, based on the total water volume input from the Sandusky River. For example, 205,367 particles are released in March due to the highest river discharge in this month, which approximately equals the total number of particles (207,050) released from April to October.

Sandusky Bay Biological Model
The biological model used in this work is a general 3-D NPZD model [17]. As a common approach, the biological model is constructed by implementing 1-D NPZD models for each vertical column of PCPM cells that are distributed spatially across the 2-D domain to form a 3-D representation of the system. An exchange of properties between adjacent water columns occurs across their shared interface through advection and dispersion. Figure 4 displays the interactions among state variables in the NPZD model.

Sandusky Bay Biological Model
The biological model used in this work is a general 3-D NPZD model [17]. As a common approach, the biological model is constructed by implementing 1-D NPZD models for each vertical column of PCPM cells that are distributed spatially across the 2-D domain to form a 3-D representation of the system. An exchange of properties between adjacent water columns occurs across their shared interface through advection and dispersion. Figure 4 displays the interactions among state variables in the NPZD model.  The governing equations for the model framework are based on Luo et al. [17], and the mathematical expressions for each term of the system of Equations (10) is presented in Appendix A. Several equations in the governing equations are modified for this study based on literature review. The light-limited, nutrient-limited, and temperature-limited functions ( ), ( ), ( ), respectively, that contribute to the ( ) are taken from Platt et al. [41] and Nicklisch et al. [42]. Also, the light attenuation functions are adjusted to Rowe et al. [18].
where , are the initial linear slope at low irradiance and the negative slope at the high irradiance that characterizes photoinhibition [43], is the maximum potential growth rate, and is the light intensity. The nutrient threshold represents the pool of nutrient that was assumed to be biologically unavailable. and are the optimal growth temperature and minimal growth temperature, respectively.
is the light attenuation coefficient that accounts for the impact of water turbidity, phytoplankton, and detritus on the light attenuation. Model parameterization is based on literature review [17,18,41,43] and subjective tuning for the Sandusky Bay simulation as there is no established NPZD model for the Sandusky Bay region. (See Table 2 for model parameterization).

Idealized Case 1: Advection-Diffusion Plume
To illustrate the effect of the cell-based averaging (step 6), we show results of the first idealized case for four different values of the cell-based redistribution parameter ( = 0, 0.01, 0.1, 0.5) in Figure  5. In Figure 5, there are three panels for each value of . The top panel shows the locations of particles after 720 time steps (12 min). The particles are colored using a blue-to-red scale for concentration values from 0 to 1. Particles with a concentration value of exactly 0 are colored light gray. The second The governing equations for the model framework are based on Luo et al. [17], and the mathematical expressions for each term of the system of Equations (10) is presented in Appendix A. Several equations in the governing equations are modified for this study based on literature review. The light-limited, nutrient-limited, and temperature-limited functions f (I), f (N), f (T), respectively, that contribute to the P(uptake) are taken from Platt et al. [41] and Nicklisch et al. [42]. Also, the light attenuation functions are adjusted to Rowe et al. [18].
where α I , β I are the initial linear slope at low irradiance and the negative slope at the high irradiance that characterizes photoinhibition [43], µ max is the maximum potential growth rate, and I is the light intensity. The nutrient threshold N 0 represents the pool of nutrient that was assumed to be biologically unavailable. T opt and T min are the optimal growth temperature and minimal growth temperature, respectively. k d is the light attenuation coefficient that accounts for the impact of water turbidity, phytoplankton, and detritus on the light attenuation. Model parameterization is based on literature review [17,18,41,43] and subjective tuning for the Sandusky Bay simulation as there is no established NPZD model for the Sandusky Bay region. (See Table 2 for model parameterization).

Idealized Case 1: Advection-Diffusion Plume
To illustrate the effect of the cell-based averaging (step 6), we show results of the first idealized case for four different values of the cell-based redistribution parameter (α = 0, 0.01, 0.1, 0.5) in Figure 5. In Figure 5 [44,45], i.e., where C(x) is the centerline concentration x meters away from the channel entrance. erf represents Gauss error function, defined as: erf(x) = 2 √ π x 0 e −t 2 dt. In the case α = 0, there is no cell-based redistribution of properties, so all particles retain their initial concentration values of either C = 0 (light gray in panel 1) or C = 1 (red in panel 1). As seen in the second and third panels, the random-walk diffusion in the particle trajectories does provide a rough approximation to the analytical solution by mixing of C = 0 and C = 1 particles in the PCPM cells. Of course increasing the number of particles in the simulation would provide a more accurate approximation, but would also increase the computational load. Setting the cell-based redistribution parameter to even the small value of α = 0.01 provides a significant improvement in the solution with the same number of particles, particularly for x > 500 m. Now particles can have any value of C between 0 and 1. Increasing the redistribution parameter to α = 0.1 further improves the solution for x < 500 m. Further increasing α to 0.5 does not significantly improve the solution in comparison to α = 0.1.

Idealized Case 2: Vertical Settling
The results at the end of 5000 time steps of the second idealized case are shown in Figure 6. In Figure 6, the dots represent the locations of the particles on the vertical axis and the value of concentration they are carrying on the horizontal axis. The thin line is the cell average concentration. The thick line is the analytical solution, C = e −ws kz z (16) As shown in Figure 6, the model properly simulates the change in concentration due to vertical settling and mixing while allowing the particles to remain approximately uniformly distributed in the vertical. The simulation accuracy increases with increased resolution of vertical layers. The model result with 20 vertical layers shows a close agreement with the analytical solution. Specifically, Figure 7 shows the evolution in time of the root mean square difference (RMSD) between the cell averages and the analytical solution for the three cases. While the RMSD in the simulation with 5 layers remains above 0.2 (the magnitude of initial error) over the entire simulation, the RMSD decreases quickly to 0.02 after 500 time steps and remains stable at this level when vertical resolution is increased to 20 layers.

Idealized Case 2: Vertical Settling
The results at the end of 5000 time steps of the second idealized case are shown in Figure 6. In Figure 6, the dots represent the locations of the particles on the vertical axis and the value of concentration they are carrying on the horizontal axis. The thin line is the cell average concentration. The thick line is the analytical solution, = As shown in Figure 6, the model properly simulates the change in concentration due to vertical settling and mixing while allowing the particles to remain approximately uniformly distributed in the vertical. The simulation accuracy increases with increased resolution of vertical layers. The model result with 20 vertical layers shows a close agreement with the analytical solution. Specifically, Figure  7 shows the evolution in time of the root mean square difference (RMSD) between the cell averages and the analytical solution for the three cases. While the RMSD in the simulation with 5 layers remains above 0.2 (the magnitude of initial error) over the entire simulation, the RMSD decreases quickly to 0.02 after 500 time steps and remains stable at this level when vertical resolution is increased to 20 layers.

Application to Sandusky Bay
To ensure the validity of the 1-D NPZD biological model, it was configured to duplicate several scenarios (not shown) from Edwards et al. [46]. The model demonstrated the expected linear stability of a vertically-distributed, NPZ ecosystem model when it was used in the scenarios that incorporate the impact of vertical mixing on biological dynamics. Scenarios include stable profiles, damping oscillatory dynamical trajectories, and vertically phase-locked systems, depending on the depth and choice of parameters and strength of vertical diffusion, which can be discerned from the eigenvalues in linear stability analysis [46]. The 1-D NPZD model used in this study reproduced all of these cases almost identically.
Before examining the impact of physical transport on the biological dynamics in Sandusky Bay, we first tested the ability of PCPM to simulate the advection-diffusion using a conservative soluble tracer in a natural setting. The Sandusky River plume is simulated using the conventional solubletracer model based on the advection-diffusion equation and compared to the result of PCPM approach in Figure 8. It is clear that the plumes simulated using the two methods show a very similar pattern, indicating the validity of the PCPM. Upon closer review, the plume simulated with solubletracer model shows a smoother evolution and sharper gradient near the plume front. We speculate this is partly due to the constant random-walk scale (10 m /s) used in the current particle-tracking model configuration. Nonetheless, the attractiveness of PCPM is its computational efficiency; it runs ~100 times faster than the soluble-tracer model which will be discussed in detail in the following section and Table 2.

Application to Sandusky Bay
To ensure the validity of the 1-D NPZD biological model, it was configured to duplicate several scenarios (not shown) from Edwards et al. [46]. The model demonstrated the expected linear stability of a vertically-distributed, NPZ ecosystem model when it was used in the scenarios that incorporate the impact of vertical mixing on biological dynamics. Scenarios include stable profiles, damping oscillatory dynamical trajectories, and vertically phase-locked systems, depending on the depth and choice of parameters and strength of vertical diffusion, which can be discerned from the eigenvalues in linear stability analysis [46]. The 1-D NPZD model used in this study reproduced all of these cases almost identically.
Before examining the impact of physical transport on the biological dynamics in Sandusky Bay, we first tested the ability of PCPM to simulate the advection-diffusion using a conservative soluble tracer in a natural setting. The Sandusky River plume is simulated using the conventional soluble-tracer model based on the advection-diffusion equation and compared to the result of PCPM approach in Figure 8. It is clear that the plumes simulated using the two methods show a very similar pattern, indicating the validity of the PCPM. Upon closer review, the plume simulated with soluble-tracer model shows a smoother evolution and sharper gradient near the plume front. We speculate this is partly due to the constant random-walk scale (10 m 2 /s) used in the current particle-tracking model configuration. Nonetheless, the attractiveness of PCPM is its computational efficiency; it runs~100 times faster than the soluble-tracer model which will be discussed in detail in the following section and Table 2. To aid in model development, several datasets are gathered from literature as well as data acquisition organizations. Sandusky River daily discharge and nitrogen concentration are available from the National Center for Water Quality Research (https://ncwqr.org/monitoring/data/, accessed on 12 July 2018). Nitrogen, Chlorophyll concentration, and in-situ temperature data are available from two observational sites (ODNR1 and EC1163) in the eastern bay from May-October 2015, sampled by Bowling Green State University [28].
Using the PCPM-NPZD model, the importance of physical transport is demonstrated by comparing model results between the PCPM-NPZD-NOADV (no advection) simulation and the realistic PCPM-NPZD-REAL simulation. In the PCPM-NPZD-NOADV simulation, the model is configured the same as the PCPM-NPZD-REAL, except the movement of particles is driven only by turbulent diffusion without including advective processes due to the Sandusky River and wind field. Each simulation consists of three cases that use a high initial nutrient concentration of 9 mg N/L in June (case 1), medium initial concentration of 0.46 mg N/L averaged from July to September (case 2), and low initial concentration of 0.0075 mg N/L in August (case 3), respectively. The concentration values are estimated as total nutrient loading from the Sandusky River to the Bay (https://ncwqr.org/monitoring/data/, accessed on 12 July 2018) divided by the total water volume of the bay. The comparison of model results is presented in Figure 9. The simulation of PCPM-NPZD-NOADV without resolving the advective transport processes shows a significant discrepancy from observational data (Figure 9, upper panels). The model fails to capture both the timing and magnitude of the blooms in all three cases, and model results are sensitive to the initial nutrient concentration.
On the other hand, after the impact of advective processes is resolved in the PCPM-NPZD-REAL simulation, the model accurately depicts the magnitude of the chlorophyll peak in mid-August (Figure 9, lower panels), and model results are insensitive to the initial condition, but determined by To aid in model development, several datasets are gathered from literature as well as data acquisition organizations. Sandusky River daily discharge and nitrogen concentration are available from the National Center for Water Quality Research (https://ncwqr.org/monitoring/data/, accessed on 12 July 2018). Nitrogen, Chlorophyll concentration, and in-situ temperature data are available from two observational sites (ODNR1 and EC1163) in the eastern bay from May-October 2015, sampled by Bowling Green State University [28].
Using the PCPM-NPZD model, the importance of physical transport is demonstrated by comparing model results between the PCPM-NPZD-NOADV (no advection) simulation and the realistic PCPM-NPZD-REAL simulation. In the PCPM-NPZD-NOADV simulation, the model is configured the same as the PCPM-NPZD-REAL, except the movement of particles is driven only by turbulent diffusion without including advective processes due to the Sandusky River and wind field. Each simulation consists of three cases that use a high initial nutrient concentration of 9 mg N/L in June (case 1), medium initial concentration of 0.46 mg N/L averaged from July to September (case 2), and low initial concentration of 0.0075 mg N/L in August (case 3), respectively. The concentration values are estimated as total nutrient loading from the Sandusky River to the Bay (https://ncwqr.org/monitoring/data/,accessedon12July2018) divided by the total water volume of the bay. The comparison of model results is presented in Figure 9. The simulation of PCPM-NPZD-NOADV without resolving the advective transport processes shows a significant discrepancy from observational data (Figure 9, upper panels). The model fails to capture both the timing and magnitude of the blooms in all three cases, and model results are sensitive to the initial nutrient concentration.
On the other hand, after the impact of advective processes is resolved in the PCPM-NPZD-REAL simulation, the model accurately depicts the magnitude of the chlorophyll peak in mid-August ( Figure 9, lower panels), and model results are insensitive to the initial condition, but determined by Sandusky River discharge and its nutrient loading. Three cases show nearly identical results. Results also support the field sampling study in Salk et al. [28]. Their study finds a strong, non-linear connection between the bloom occurrence and the hydraulic residence time of the bay, which varies dramatically from 8 days to several months depending on the Sandusky River flow rate in the physical transport process. Although further development of the NPZD is undoubtedly necessary to resolve the onset and variability of the algal blooms by refining the structure of the biological model and improving the model parameterization, it is beyond the scope of this work in which we focus on demonstrating the feasibility of linking hydrodynamic effects and biological processes through the PCPM in a Lagrangian/Eulerian framework. Further development of the biological model and its application to the study of mechanisms responsible for the HABs in Sandusky Bay will be presented in another paper. Sandusky River discharge and its nutrient loading. Three cases show nearly identical results. Results also support the field sampling study in Salk et al. [28]. Their study finds a strong, non-linear connection between the bloom occurrence and the hydraulic residence time of the bay, which varies dramatically from 8 days to several months depending on the Sandusky River flow rate in the physical transport process. Although further development of the NPZD is undoubtedly necessary to resolve the onset and variability of the algal blooms by refining the structure of the biological model and improving the model parameterization, it is beyond the scope of this work in which we focus on demonstrating the feasibility of linking hydrodynamic effects and biological processes through the PCPM in a Lagrangian/Eulerian framework. Further development of the biological model and its application to the study of mechanisms responsible for the HABs in Sandusky Bay will be presented in another paper.

Summary and Conclusions
In this paper, we describe a novel method by integrating a property-carrying particle model (PCPM) and an Eulerian concentration biological model for ecosystem modeling. The model is tested in idealized cases and its utility is demonstrated in a practical application to Sandusky Bay. The novelty of this new technique lies in its integration of hydrodynamic effects via the property-carrying particle tracking model and Eulerian grid-based biological modeling approach. Overall, there are several advantages of the PCPM over traditional Eulerian-based tracer approaches. The PCPM is simpler to implement and more efficient as it does not need to solve the advection-diffusion equation. Instead, the PCPM uses pre-computed particle trajectories to resolve the hydrodynamic condition based on currents from a hydrodynamic model. This means that the hydrodynamic model only needs to be run once giving one the ability to run different biological scenarios for the same physical characteristics; ultimately saving significant computational time.
As summarized in Table 2, a 30-day hydrodynamic simulation for the Lake Erie-Sandusky Bay FVCOM model takes 15 h to complete using 64 central processing units (CPUs). This step is necessary for both the PCPM-NPZD model and the traditional Eulerian, grid-based biophysical model. Once the hydrodynamic simulation is complete, PCPM-NPZD needs to calculate the Lagrangian trajectories of a large number of current-following tracer particles; it takes 1.5 h using a single CPU to track~290,000 particles in the domain in March, which tracks the largest number of particles in the bay within a single month. The PCPM-NPZD model simulation using the particle trajectories as input completes a 30-day simulation within 5 min using a single CPU while it takes 5 h for a traditional Eulerian, grid-based biophysical model to complete the same simulation using 32 CPUs. If compared with the same computational power (e.g., single CPU), the PCPM-NPZD approach runs~100 times faster for the biophysical modeling ( Table 2; Scenario 1).
More importantly, in the PCPM framework, the hydrodynamics and associated water transport and mixing represented by particle trajectories are "reserved" and not affected by biochemical properties. In other words, it only takes another 5 min to run the PCPM-NPZD for a different set of parameters and property configurations. This is extremely useful during the model calibration and or ensemble simulations. The PCPM-NPZD would take 11.5 h to complete 100 runs using a single CPU, while the traditional method would require 500 h simulation using 32 CPUs (Table 2; Scenario 2). Such a high level of efficiency is not available from tracer-based models because one will have to re-run the Eulerian-based biological model for any change in parameter configuration or estimation of different property concentration. Also, the PCPM is capable of providing comparable simulation results to the soluble-tracer model, although the global and local mass conservation is not strictly preserved with finite particles. Above all, it is the PCPM's computational efficiency and coupling flexibility which makes it an attractive alternative method to the traditional approach. Author Contributions: Conceptualization, P.X. and D.J.S.; Formal analysis, X.Z., C.H., R.K. and X.Y.; Project administration, P.X.; Validation, X.Z.; Visualization, X.Z., C.H. and X.Y.; Writing-original draft, P.X. and D.J.S.; Writing-review and editing, P.X., D.J.S., R.K.
Funding: Xue's work is partly supported by the Ohio Department of Natural Resources, Sandusky Bay Initiative.
Acknowledgments: This is the contribution 57 of the Great Lakes Research Center at Michigan Technological University. The Michigan Tech high performance computing cluster, Superior, was used in obtaining the modeling results presented in this publication. We thank George Bullerjahn, Robert McKay and Timothy Davis (Bowling Green State University) for access to their Sandusky Bay water quality data.

Conflicts of Interest:
The authors declare no conflict of interest. The mathematical expressions for the biological terms in the system of Equation (10) are listed below, where the definition and value of each parameter in the equations is above.