Numerical Simulation on the Diffusion of Alien Phytoplankton in Bohai Bay

In order to investigate the motion feature of alien phytoplankton in the Bohai Bay area, this manuscript builds a two-dimensional tide hydrodynamic model coupled with a particle tracking model to simulate the alien phytoplankton movement trajectory and the diffusion processes with different specific growth rates in three major ports of Bohai Bay (Caofeidian port, Tianjin port, and Huanghua port). The results show that the movement of alien phytoplankton is mainly affected by the tidal circulation near the port, and the diffusion trend is basically consistent with the residual flow in Bohai Bay. The distribution density of alien phytoplankton is directly affected by the specific growth rate of the population and is positively related to specific growth rate. The released alien phytoplankton in the three major ports are all concentrated around the ports area. The largest distribution density is in Tianjin port, and the possibility of red tide disasters is also greatest here compared with the other two major ports. It is necessary to strengthen the monitoring of alien organisms in the port area and actively prevent alien phytoplankton from entering Bohai Bay through ship ballast water.


INTRODUCTION
Bohai Bay is an important economic belt in China, and it is located in the west of Bohai Sea, with unique geographical and important strategic position. There are many ports in Bohai Bay, and some non-inactivated alien phytoplankton (sporangia) may be carried in with ships arriving at the port. Since 2012, Bohai Bay has been the most seriously affected area by foreign phytoplankton (State Oceanic Administration, 2019). Up to 2019, 17 species of alien phytoplankton have been found in Bohai Bay, among which 16 species belong to red tide organisms (Pu et al., 2020). Under the action of suitable hydrometeorology, seawater physical and chemical factors, water nutrients, and other factors, red tide will be triggered, which will cause a serious negative impact on the marine environment and social economy.
Many scholars studied the ecological environment of Bohai Bay, and phytoplankton has become the focus of research. In July 2003, Cao et al. (2004Cao et al. ( , 2006 detected 31 species of phytoplankton belonging to 19 genera and found red tide organisms in five sites near Luju River in Tianjin port. Wu et al. (2006), by constructing a three-dimensional dynamic ecological coupling model, found that tidal current has an impact on the vertical distribution of phytoplankton and can promote the formation of red tide when the light and nutrient conditions are appropriate, which provides a dynamic basis for further research on the occurrence of red tide. Liu (2007) summarized the temporal and spatial distribution characteristics of phytoplankton quantity and the distribution of dominant species by analyzing phytoplankton data in Bohai Bay in 2005. The number of phytoplankton cells, the diversity of phytoplankton, and the complexity of community structure in Bohai Bay showed seasonal variation Yang et al., 2007). Yin et al. (2007) combined the survey results of Tianjin coastal zone and the ecological survey data of Tianjin coastal area in Bohai Bay from 2008 to 2012 and determined the species composition of phytoplankton in summer in Tianjin coastal area of Bohai Bay. In addition, many studies have shown that the community structure of phytoplankton is affected by multiple environmental factors (Liu X. T., 2011;Yin et al., 2013aYin et al., ,b, 2015Wang, 2015;Wang et al., 2016;Zhang et al., 2016). Qin et al. (2017) proposed that there are five kinds of alien phytoplankton in Bohai Bay through a large number of research and investigation, and evaluated their ecological risk by constructing an evaluation model. Liu (2019) studied the temporal and spatial distribution characteristics and response relationship of nutrients and phytoplankton communities in Bohai Bay, which provided a scientific basis for further understanding the evolution characteristics of the ecological environment in Bohai Bay under the joint action of human activities and climate change. Pang et al. (2020) constructed the eutrophication evaluation index system of Tianjin coastal waters based on the PSR model, evaluated and verified the eutrophication level of Tianjin coastal waters by using the analytic hierarchy process, and explained that human behavior had a negative impact on the ecosystem of Tianjin coastal waters in Bohai Bay.
At present, most studies on alien phytoplankton in the world stay on their own biological characteristics (Goes et al., 2013;Cristian et al., 2016;Silkin et al., 2016;Belyaeva, 2019), and there are few studies on their transport and diffusion under the influence of hydrodynamic conditions; the relevant studies on phytoplankton in Bohai Bay are also based on the research of Chinese scholars.
The application of numerical simulation technology to predict the movement and diffusion of invasive phytoplankton in Bohai Bay is of great significance in the monitoring and prevention of red tide disaster in Bohai Bay. The particle tracking model is based on the hydrodynamic model and superposed with the particle tracking module to study the transport and diffusion process of particles under the influence of wind and hydrodynamic conditions. The particle tracking model is very suitable for establishing the tracking model and simulating the particle track and tracking trajectory in the water environment. Pablo (2005) used the Lagrange particle tracking method to simulate the movement characteristics of eel population in the upper reaches of the river. Zhang et al. (2012) used the Lagrange method to track the path of jellyfish based on the simulation results of the sea flow field. Han et al. (2013) based on the results of a three-dimensional hydrodynamic model of a puppet lake, using the particle tracking method to simulate the movement characteristics of algae, ignoring the growth and reproduction of algae, and studied the accumulation of algae under different wind directions. Zheng (2014) simulated the drift path and diffusion range of planktonic larvae under the dual influence of wind field and tidal current field based on the three-dimensional hydrodynamic model of the Yellow Sea area and superimposed particle tracking module, and the simulated data were basically consistent with the measured data. The particle tracking method is widely used to study the movement characteristics and migration path of alien phytoplankton, and the results are reasonable and satisfactory. However, previous studies on plankton simulation ignored the growth and reproduction of plankton.
In view of the above shortcomings, based on the numerical simulation of tidal current field in a specific period of Bohai Bay, combined with the residual current theory and particle tracking model, the tide-induced Euler residual current field and tide-induced Lagrangian residual current field are analyzed and compared; combined with the biological characteristics of alien phytoplankton, the Lagrangian particle tracking model is applied to analyze the distribution of alien phytoplankton (sporangia) in three ports of Bohai Bay in order to provide support for the prevention and control of red tide disaster in Bohai Bay.

Numerical Simulation of the Tidal Current Field in Bohai Bay
In order to conduct a numerical simulation of the tidal current field in Bohai Bay (the area surrounded by CD open boundary and land boundary in Figure 1), this study intends to build a tidal hydrodynamic numerical model of Bohai Sea (the area surrounded by AB open boundary and land boundary in Figure 1) with a wider coverage and then provide water level for the CD open boundary required by the numerical model of Bohai Bay.
The period of hydrodynamic simulation in Bohai Sea is from April 15 to May 15, 2014, and the simulation time step is 200 s. The model is calibrated by using the data of tidal table water level at four tidal stations in Bohai Sea, and the current field in Bohai Sea can be properly simulated (Pu, 2019) when the eddy viscosity coefficient C s = 0.28, Manning coefficient M = 48 m 1/3 /s, i.e., n = 0.0208.
The current model encrypts the calculation grid of the offshore area, with 6,316 nodes and 10,335 grids in the whole Bohai Bay simulation area. The grid division is shown in Figure 2A. The open boundary tidal level data of Bohai Bay are obtained from the results of the Bohai Sea current hydrodynamic model mentioned above, and the simulation period and parameter settings are the same as those of the Bohai Sea current model. The simulation results are verified by the measured data of tide level and tidal current stations in Bohai Bay. As shown in Figure 2B, T1 and T6 are tide level stations and S1 and S3 are tidal current stations. The simulated results were verified by using the observed values of sea level, velocity, and flow direction of four stations from April 30, 2014, 10:00 to May 1, 2014 10:00 (spring tide period) and from May 12, 2014, 07:00 to May 13, 2014 07:00 (neap tide period).  The results of Figures 3, 4 show that the simulated current field in Bohai Bay is reasonable and feasible, and further simulation work can be carried out based on the flow field results.

Analysis and Comparison of Tide-Induced Residual Current
The non-linear effect of tidal current will produce tide-induced residual current, which will further form tidal circulation and plays an important role in the migration and distribution of material in the sea area. The tide-induced residual current can be described by two methods: Euler residual current and Lagrangian residual current. Euler residual current is the average value of instantaneous velocity at fixed position in the flow field in certain tidal cycles; Lagrangian residual current is the ratio of the displacement of fluid microcluster to the time consumed in a certain tidal cycle .

Analysis of Euler Residual Current Field
According to the definition of Euler residual current and based on the hydrodynamic field simulation results above, the  tide-induced Euler residual current in Bohai Bay is obtained as shown in Figure 5. The chart shows that the overall trend of the residual current in Bohai Bay is that seawater flows into Bohai Bay along the coast and flows out from the middle of the open boundary to maintain water balance. The residual current in the middle of Bohai Bay is weak, and it is strong along the coast. The residual current near the coast of Caofeidian port, Dagukou, Dakou River, and Dongfeng port is obviously enhanced. A clockwise circulation is between Dagukou and the west of Caofeidian port, and a counterclockwise circulation is in the southwest of Bohai Bay. A counterclockwise circulation forms in Caofeidian port, and in the south of Caofeidian port another clockwise circulation forms, and the two circulations present a double circulation structure.

Analysis of Lagrangian Residual Current Field
The Lagrangian particle tracking model is constructed to simulate the characteristics of Lagrangian residual current in Bohai Bay. A certain number of tracer particles are released at three locations as shown in Figure 6. Under the same hydrodynamic conditions, the particle migration path trajectory is shown in Figure 7. The simulation results show that after the particles are released in Caofeidian port, they move mainly in the east-west direction with the longshore current and do not expand to the central  area of Bohai Bay. In Tianjin port, the particles released mainly diffuse around the port in the simulation time, and the motion track basically covers the area of Tianjin port. After being released in Huanghua port, the particles move back and forth with the longshore current, and they mainly diffuse along the coast and move to the northwest of Bohai Bay.
Theoretically, with first-order approximation, the Lagrangian residual current is the material transport velocity, and the difference between the Euler residual current and the material transport velocity is a Stokes drift velocity; therefore, the Euler residual current and the Lagrangian residual current are not completely consistent. In practical application, it is generally considered that the difference between the two is small enough to be ignored, but relevant studies show that the difference cannot be ignored under certain conditions (Wei et al., 2003), and it is necessary to use the Lagrangian method to study the material migration in the sea area (Yu and Chen, 1983;Tang, 1987). In view of this, although the shape of the Eulerian residual current field has great reference significance for predicting the motion trend and trajectory range of Lagrangian particles, the final trajectory coverage of Lagrangian particles in the three ports is highly consistent with the shape of Eulerian residual flow field in the three ports. The following work will be conducted based on the Lagrangian method to study the transport characteristics of alien organisms.

Numerical Simulation of Alien Phytoplankton Diffusion
In this study, we continue to use the Lagrangian particle tracking model to characterize the movement of alien phytoplankton by Lagrangian particle trajectory, and simulate the movement of alien phytoplankton (sporangium) in Caofeidian port, Tianjin port, and Huanghua port of Bohai Bay, revealing the movement characteristics of alien phytoplankton in Bohai Bay, so as to provide support for the prevention and control of red tide disaster in Bohai Bay.

Biological Characteristics of Alien Phytoplankton
In order to carry out the relevant simulation work, it is necessary to make clear the mathematical description of the growth, reproduction, and attenuation process of alien phytoplankton in Bohai Bay, so as to reflect its quantitative change process in water. The number of phytoplankton community in water is directly affected by the maximum specific growth rate, zooplankton feeding, competition, population migration, and other factors.

Specific growth rate of phytoplankton
The life process of phytoplankton includes growth and death, and its instantaneous number is the difference between the increment of plant cells and the decrement of plant cells (including sedimentation and death of phytoplankton) in the instantaneous interval. The net growth rate can be expressed by the specific growth rate (µ) of the number of phytoplankton, and µ can be expressed as Equation 1.
In the formula, µ (day −1 ) is the specific growth rate of phytoplankton, N 1 (cell L −1 ) is the number of phytoplankton at time t 1 , and N 0 (cell L −1 ) is the number of phytoplankton at t 0 .
µ is a non-fixed parameter, which is not only determined by the genetic material of individual cells but also affected by nutritional level, light, temperature, and other external environmental factors. When the phytoplankton population grows in the most suitable environmental conditions, µ is a constant value, called the maximum specific growth rate µ max or intrinsic growth rate (Sun and Ning, 2005). Moreover, different phytoplankton have different maximum specific growth rates.
Through the investigation of a large number of literature (Wang et al., 2001;Yu, 2005;Li, 2006;Hu et al., 2010;Liu S. X., 2011;Ma et al., 2012;Ge et al., 2017), the maximum specific growth rates of 17 kinds of alien phytoplankton in Bohai Bay were obtained. The largest intrinsic growth rate of these species was found in Pseudo-nitzschia pungens (0.92 day −1 ), and the smallest was found in Cochlodinium polykrikoides (0.42 day −1 ), with an average of 0.52 day − 1 .

Feeding rates of zooplankton
In natural waters, zooplankton mainly feed on various phytoplankton and suspended organic matter, and the feeding rate of zooplankton varies in different sea areas. Zhang and Wang (2000) sampled and analyzed the spatial distribution of microzooplankton in the Bohai Sea and calculated the feeding rate of zooplankton as 0.42-0.69 day −1 by the dilution method. In order to determine the maximum influence range of alien phytoplankton diffusion as large as possible, the minimum feeding rate of zooplankton 0.42 day −1 was taken in this paper.

Specific growth rate of the phytoplankton community
The specific growth rate of the phytoplankton community is not only directly affected by the maximum specific growth rate but also controlled by a variety of factors, including zooplankton feeding, competition, and population migration. According to the research purpose, this manuscript only took the specific growth rate of a single phytoplankton community into consideration, and factors such as population competition and population migration were not considered.
The specific growth rate of phytoplankton community can be expressed as Equation 2.
In Equation 2, µ is the community specific growth rate, µ max is the maximum specific growth rate, and g is the zooplankton feeding rate.
For 17 species of phytoplankton in Bohai Bay, the maximum specific growth rate of the community was 0.5 day −1 , the minimum was 0 day −1 , and the average was 0.1 day −1 .
Under ideal conditions, without considering the maximum environmental capacity, the growth of a single phytoplankton community conforms to the exponential growth model (Fan et al., 2011)

expressed as Equation 3.
In the formula, N t is the biomass of phytoplankton at time t, N 0 is the biomass of phytoplankton at the initial time, t is the growth time, and µ is the community specific growth rate.

Influence of temperature factor
Since the viability of phytoplankton sporangium is much better than that of phytoplankton cells, and the survival ability in ballast water is stronger, the alien phytoplankton leaked at the port is mainly the non-inactivated sporangium in ballast water. The sporangium has a certain dormant period, and it will germinate and grow as long as the environment is suitable.
The Bohai Bay region is a temperate continental climate, with low water temperature from December to March. From April to spring, the water temperature gradually increases, and the surface water temperature is high (11-12 • C) along the coast, and the offshore water temperature is low (7-8 • C). The spores of alien phytoplankton begin to germinate in this season. From April to May, with water temperature gradually increasing, it will be more suitable for the growth of alien phytoplankton, so during the whole simulation period, the alien phytoplankton is in the growth and reproduction period. Therefore, the model particles will be released throughout the simulation time.

Diffusion Simulation of Alien Phytoplankton Particlization hypothesis of phytoplankton
Considering that the alien phytoplankton in the study area are single-cell organisms and have small volume and mass, under the action of gravity, buoyancy, inertia, and other forces in the sea, they move with the tide. We can learn from the sediment movement research methods to generally analyze the movement of alien phytoplankton.
In order to explore the transport and diffusion law of alien phytoplankton under the tidal current field in Bohai Bay, the research model is appropriately generalized, and the rules are as follows: (1) Alien phytoplankton is represented by particles, and the movement of phytoplankton is ignored.
(2) Based on the two-dimensional hydrodynamic simulation results, this manuscript only considers the transportation and dispersion of particles in the horizontal direction without considering settlement.

The setting of phytoplankton growth
The maximum distance between the open boundary of Bohai Bay and the port on the land boundary is no more than six nautical miles, with an average water depth of 12.5 m. According to the relevant discharge requirements of the ballast water convention (Fu et al., 2016), under normal circumstances, the possibility of alien phytoplankton spreading from the ballast water replacement site to the Bohai Bay is small. However, in the process of ballast water treatment, due to failure of the treatment technology to inactivate all organisms (Sun, 2014), the ballast water may contain some algae cysts that have not been inactivated, and these cysts will enter the port with ships and be released at the port. Based on the research results on the invasion pathway of alien phytoplankton in Bohai Bay, the bay is an important economic zone in China, with many coastal ports, among which Caofeidian port, Tianjin port, and Huanghua port are the three largest ports in the bay, and the possibility of alien phytoplankton cyst leakage in the three ports is great. In this study, Caofeidian port, Tianjin port, and the inner bay port of Huanghua port in Bohai Bay are selected as the key research areas. Considering the large area of the port, the point source particles near the central point of the port are selected in this model. The particle release sources are shown in Figure 6.
In the model, the instantaneous particle number is the instantaneous total number of alien phytoplankton. The number of phytoplankton sporangium in the cabin is very large, which can reach hundreds of millions (Long et al., 2005). To see the results more clearly, in this study, the number of sporangium leaked at a port is assumed as 10 6 .
The 17 species of alien phytoplankton in Bohai Bay are all single-cell organisms with small size, and the mass of a single phytoplankton is about 0.1 µg, so the initial total mass released in the simulation study is 10 5 µg. The growth rates of the 17 alien phytoplankton population are 0-0.5 day −1 . In order to more accurately learn the possible diffusion range of alien phytoplankton in Bohai Bay, the transport and diffusion simulations of alien phytoplankton are taken with different growth rates of species community, which are minimum value (µ = 0 day −1 ), average value (µ = 0.1 day −1 ), and maximum value (µ = 0.5 day −1 ), respectively.
As mentioned above, the number of phytoplankton increased exponentially with time during the simulation period. Accordingly, the time series file of particle release rate changing with time should be established in the model. Figures 8-10 show changes of total mass of alien phytoplankton corresponding to different specific growth rates of phytoplankton communities in the simulation period. The mass of particles released by the particle tracking model over time should meet this rule.

The setting of simulation conditions
Based on the flow field simulation of the study area, the particle tracking model was coupled to simulate the transport and diffusion of particles released at specific locations during April 15 to May 15, 2014. At the initial time, the tracer particles are released at Caofeidian port, Tianjin port, and Huanghua port, respectively. The instantaneous trajectory was output once per hour, and 720 movement moments are recorded. The description of simulated working conditions is shown in Table 1. Figures 11-13 show the density profile of released alien phytoplankton under different working conditions in three ports after 30 days of simulation.

RESULTS
Through the simulation and analysis of particle diffusion trajectory, particle number, diffusion range, and diffusion trend in Caofeidian port, Tianjin port, and Huanghua port, the following information can be obtained: the movement of alien phytoplankton is mainly affected by the tidal circulation in the port area; for different release quantities, the diffusion trend and diffusion range of alien phytoplankton in the same area are basically the same, that is, the number of alien phytoplankton will not have a significant impact on the diffusion range in the study area, but the release quantity will have a direct impact on TABLE 1 | Working conditions setting of phytoplankton community specific growth rate in each port.

Working conditions Instruction
Condition 1 Minimum specific growth rate of alien phytoplankton community µ = 0 day −1

Condition 2
Average specific growth rate of alien phytoplankton community µ = 0.1 day −1

Condition 3
Maximum specific growth rate of alien phytoplankton community µ = 0.5 day −1 the density distribution of alien phytoplankton in the study area, and the larger the amount of phytoplankton released, the greater the phytoplankton density in the water body.

DISCUSSION
The existing application of the particle tracking model is basically based on the application background, which simplifies the research object to particles with no life characteristics (even the simulation research on the migration of aquatic community is the same), that is, the research object can only drift with the current and has no active ability. In this study, based on the investigation of the species and growth characteristics of alien phytoplankton in Bohai Bay, the concept of specific growth rate is used to endow the particle (alien phytoplankton) with "reproductive capacity, " that is, measures are taken to make the particle (alien phytoplankton) propagate with time, and its number in the water area continues to increase. The improved application of the model makes it more suitable for the actual scene under specific conditions than the conventional particle tracking model. It should be pointed out that the specific growth rate data used in the study are not measured in Bohai Bay, and there are some differences in the specific growth rate values of different researchers for the same phytoplankton. In order to ensure that the research results have a sufficient margin of safety, the specific growth rate data with larger values are selected in this study. It is necessary to measure the specific growth rate of alien phytoplankton one by one in Bohai Bay to provide basic support for more accurate research.
Considering that the scale of movement of phytoplankton is negligible compared with the tidal current field, the movement characteristics of phytoplankton such as its own movement and sedimentation are ignored in this study, and only the drift and diffusion in the horizontal direction are considered. In fact, this is also the limitation of the conventional particle tracking model, which is difficult to fully reflect the autonomous activity of aquatic organisms, and then more accurately reflect the relationship between them and the water ecology and hydrodynamic environment. To solve this problem, the model based on the theory of ecological agent (such as the agentbased model, ABM, which is developed by DHI) can be used to define intelligent individuals (such as planktonic larvae, fish, and mammals) and their states (weight, morphology, and habits) and describe the interaction and response relationship with hydrodynamic and aquatic ecological conditions. However, the application of this model needs to be based on the comprehensive and in-depth study of aquatic biological characteristics.
Some red tides contain biotoxins. If a red tide disaster occurs, on the one hand, it will endanger the life of shellfish, fish, shrimp, crab, and other marine life in Bohai Bay; disrupt the ecological balance; and cause serious economic losses to marine fisheries. On the other hand, it will also threaten human health and life safety. If aquatic organisms containing red tide biotoxins are eaten, they may cause poisoning or even death. Thus, it is necessary to simulate and predict the movement trajectory of invasive phytoplankton and the areas where red tides may occur   and then prevent and curb the development of red tide organisms as soon as possible. These measures are of great significance to public health management.

CONCLUSION
Under the condition of neglecting the movement of phytoplankton and sedimentation, tracer particles were used to represent alien phytoplankton in this study. Based on the hydrodynamic numerical simulation technology of Bohai Bay, the movement law of alien phytoplankton in the three ports of Bohai Bay was explored through particle transport and diffusion simulation under various conditions. The conclusions and suggestions are as follows: (1) The movement of phytoplankton in Bohai Bay is mainly affected by the hydrodynamic condition of the region, namely, tidal circulation. Influenced by the hydrodynamic force in Bohai Bay, the alien phytoplankton will be transported and spread around the port after entering Caofeidian port, Tianjin port, and Huanghua port. The floating plants in Caofeidian port have the largest influence range, and the influence range of Tianjin port is the smallest, but more is still gathered around the core area of each port.
(2) The specific growth rate of phytoplankton will directly affect the distribution density of alien phytoplankton entering Caofeidian port, Tianjin port, and Huanghua port. The higher the specific growth rate, the higher the density of alien phytoplankton in the water. (3) In view of the fact that it is difficult for the alien phytoplankton to leave the Bohai Bay area driven by hydrodynamic conditions after entering the three port areas, and with the increase of the number of reproduction, and under the effect of the external environment (nutrient content in the water, hydrometeorological conditions, seawater physical and chemical factors, etc.), the possibility of subsequent red tide disasters increases correspondingly. It is necessary to strengthen the monitoring of alien organisms in the port area and actively prevent alien phytoplankton from entering Bohai Bay through ship ballast water. (4) The specific growth rates of 17 alien phytoplankton species (including 16 species of red tide organisms) in Bohai Bay may be measured through experiments, and the red tide occurrence areas of each species in this area can be simulated and predicted by combining with the numerical model in Bohai Bay, so as to provide more targeted support for the monitoring, prevention, and ballast water treatment of red tide disasters in this area and better provide reference for public health management services for the society.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
BZ wrote the manuscript. AP provided review and comments. BZ and CX analyzed the data. QW was responsible for the methods used in the study. AP and PJ conceived the research idea of this study and provided review and comments. WT revised the content of the article and was responsible for the overall submission. All the authors were engaged in the final manuscript preparation and agreed to the publication of this manuscript.