Mechanisms of Cross-Shore Transport and Spatial Variability of Phytoplankton on a Rip-Channeled Beach

We investigated whether cross-shore distributions of coastal phytoplankton to the surf zone are controlled by hydrodynamics and their biological characteristics. Data from a rip-channeled beach indicate that concentrations of phytoplankton are higher in the surf zone than offshore. To examine how phytoplankton is transported toward the shore, we used a coupled biophysical model, comprised of a 3D physical model of coastal dynamics and an individual-based model (IBM) for tracking phytoplankton on the rip-channeled beach. Waves and wind in the biophysical model were parameterized by the conditions during the sampling period. Previous studies indicated that growth rates of phytoplankton can be enhanced by high turbulence, which might contribute to high phytoplankton concentration in the surf zone. Some numerical and laboratory works showed that turbulence can also increase the downward velocity of phytoplankton, which could be carried by onshore bottom currents and remain in the surf zone. Furthermore, we adapted the IBM with the theoretical model of diurnal vertical migration (DVM) for phytoplankton. The theoretical DVM works as follows: in the morning, phytoplankton cells adhere to air bubbles and stay at the surface and close to the shore in the daytime because onshore wind and surface current direction is usually onshore; in the late afternoon, the cells switch their attachment from air bubbles to sand grains and sink to the bottom where the water ﬂow is normally onshore at night. Finally, depth-varying growth of phytoplankton was also incorporated into the DVM module. Simulations using neutral passive particles do not give the expected results of observed patterns. All tested mechanisms, i.e., wind- and wave-driven currents, rip-current circulation, turbulence-driven growth and sinking, DVM, and depth-varying growth, enhanced onshore phytoplankton migration and cell concentrations in the surf zone, indicating that both biological traits and physical factors can be essential to phytoplankton cross-shore transport and spatial variability. Our model is open to be modiﬁed and re-parameterized, followed by further analysis and validation, so that it can be more adequate for ecological assessment of coastal areas.

We investigated whether cross-shore distributions of coastal phytoplankton to the surf zone are controlled by hydrodynamics and their biological characteristics. Data from a rip-channeled beach indicate that concentrations of phytoplankton are higher in the surf zone than offshore. To examine how phytoplankton is transported toward the shore, we used a coupled biophysical model, comprised of a 3D physical model of coastal dynamics and an individual-based model (IBM) for tracking phytoplankton on the rip-channeled beach. Waves and wind in the biophysical model were parameterized by the conditions during the sampling period. Previous studies indicated that growth rates of phytoplankton can be enhanced by high turbulence, which might contribute to high phytoplankton concentration in the surf zone. Some numerical and laboratory works showed that turbulence can also increase the downward velocity of phytoplankton, which could be carried by onshore bottom currents and remain in the surf zone. Furthermore, we adapted the IBM with the theoretical model of diurnal vertical migration (DVM) for phytoplankton. The theoretical DVM works as follows: in the morning, phytoplankton cells adhere to air bubbles and stay at the surface and close to the shore in the daytime because onshore wind and surface current direction is usually onshore; in the late afternoon, the cells switch their attachment from air bubbles to sand grains and sink to the bottom where the water flow is normally onshore at night. Finally, depth-varying growth of phytoplankton was also incorporated into the DVM module. Simulations using neutral passive particles do not give the expected results of observed patterns. All tested mechanisms, i.e., wind-and wave-driven currents, rip-current circulation, turbulence-driven growth and sinking, DVM, and depth-varying growth, enhanced onshore phytoplankton migration and cell concentrations in the surf zone, indicating that both biological traits and physical factors can be essential to phytoplankton cross-shore transport and spatial variability. Our model is open to be modified and re-parameterized, followed by further analysis and validation, so that it can be more adequate for ecological assessment of coastal areas.

INTRODUCTION
Phytoplankton dynamics in coastal water largely influence marine ecosystems, fisheries, and coastal communities. For example, with favorable environmental conditions, phytoplankton can overgrow and cause harmful algal blooms (HABs) that often increase levels of toxic substances in the coastal area and the toxins accumulate in intertidal animals, and humans can suffer poisoning by consuming those toxic animals (Landsberg, 2002). Since many HAB species reside offshore, it is important to understand how they are transported to shore and possibly enter the surf zone, which is considered a "semi-permeable barrier" (Rilov et al., 2008;Shanks et al., 2010). As Shanks et al. (2016) showed, surf zone hydrodynamics affect onshore transport of Pseudo-nitzschia, one of the species causing HABs. Recent studies revealed that onshore larval migration can be influenced by surf zone hydrodynamics and characteristics of larvae (Fujimura et al., 2014(Fujimura et al., , 2017Shanks et al., 2015Shanks et al., , 2017Shanks et al., , 2018Morgan et al., 2016Morgan et al., , 2017; however, mechanisms of onshore transport of phytoplankton are not well understood. Phytoplankton cannot be modeled as simple passive particles; some phytoplankton species can float, sink, or swim (Smayda, 1970(Smayda, , 2010. In turbulent conditions, larvae of various intertidal invertebrate species sink to the bottom (Denny and Shibata, 1989;Butman, 1990;Fuchs et al., 2004;Roy et al., 2012) or actively move downwards (Fuchs et al., 2013) that can be essential traits for onshore larval migration as a model study showed (Fujimura et al., 2014). Likewise, sinking velocity of phytoplankton cells may be increased by turbulence possibly with a different mechanism (Ruiz et al., 2004;Macías et al., 2013). This might be due to preferential downward movement of particles along the peripheries of local vortical structures (Wang and Maxey, 1993). With the same vortical mechanism, the floating velocity of positively buoyant phytoplankton cells may be enhanced by turbulence. Additionally, some phytoplankton species swim with flagella or cilia. These floating and swimming behaviors were not considered in our study.
Spatial distributions of phytoplankton in the coastal area may not be controlled only by transport, but also cell growth. Variability in growth rates of phytoplankton under different environments can make their distribution patterns more complicated. For instance, Savidge (1981) showed high turbulence could increase phytoplankton growth rate about 25-75%. This may have resulted from higher levels of exposure to light and/or nutrient uptake in turbulence.
There is also diurnal periodicity in phytoplankton. Talbot et al. (1990) proposed the following diurnal vertical migration (DVM) model of surf zone phytoplankton species based on the series of their observational studies. At night, phytoplankton cells stay on the bottom by attaching to sediment with mucus they produce. In the morning, they lose the mucus and enter the water column by turbulent flows in the surf zone. Within the water column of the surf zone, they adhere to air bubbles produced by breaking waves, and increase their buoyancy. The cells float on the surface and grow in the daytime. During late afternoon, the cells start developing the mucus again and attach to sand grains to sink to the bottom. Mechanisms of formation and shedding of mucus as well as cells attaching to and detaching from air bubbles or sand grains have not been well studied. An ability to float at the surface by attaching to air bubbles or forming a semi-stable foam consisting of many small bubbles can be found in surf zone diatom species (Lewin and Schaefer, 1983). We combined the DVM model with a typical diurnal land-sea breeze cycle (i.e., onshore wind during the day and offshore wind at night). We hypothesized that phytoplankton cells stay near the surface and are kept in the surf zone by the sea breeze-driven currents during the daytime, and stay on the bottom where they are pushed toward shore by wave-driven bottom currents in the nighttime (Figure 1). Wave stress in the bottom boundary layer induces so-called benthic streaming, a bottom current flowing in the direction of wave propagation (Longuet-Higgins, 1953). Plankton near the bottom may be transported shoreward by streaming (Fujimura et al., 2014). Streaming may be weakened during daytime onshore wind events in order to balance with onshore surface flow (Figure 1).
Beach morphology also affects cross-shore transport of plankton (Shanks et al., 2010(Shanks et al., , 2018Fujimura et al., 2013;Morgan et al., 2017). Beaches can be classified as a dissipative beach with a wide surf zone and flat slope, intermediate, or reflective beach with a narrow surf zone and steep slope (Wright and Short, 1984;McLachlan and Brown, 2006). Surf zone phytoplankton taxa are found only at dissipative and intermediate beaches, but not in more reflective surf zones (Lewin and Schaefer, 1983;Talbot et al., 1990;Shanks et al., 2018). In this study, we focused on an intermediate beach where cross-shore exchange is enhanced by rip currents. There are several types of rip currents (Dalrymple et al., 2011;Castelle et al., 2016), but here we considered bathymetrically-controlled rip currents (Figure 2). When waves approach perpendicular to shore, onshore currents converge toward shoals due to wave refraction, then alongshore feeder currents accumulated in rip channels, causing offshore-directed rip currents. Rip currents weaken offshore and merge with onshore currents. This pattern often forms rip circulations that entrain and concentrate phytoplankton cells (Talbot and Bate, 1987).
A field study at a rip-channeled beach showed that concentrations of coastal phytoplankton taxa are much higher in the rip channels than offshore or over the shoals (Shanks et al., 2018). It is very likely that surf zone hydrodynamics influence such spatial variability; however, characteristics of phytoplankton should also affect the concentration of phytoplankton within the surf zone, since completely passive particles cannot transport toward shore easily (Fujimura et al., 2014).
This study was to uncover the role of physical forcing and ecological traits in coastal phytoplankton dynamics at a relatively small scale. Here, we added theoretical parameters to our model of cross-shore phytoplankton transport. We tested the effects of turbulence-driven sinking and growth within the surf zone on onshore phytoplankton transport and the concentration of phytoplankton within the surf zone. Our model built upon that proposed by Talbot et al. (1990) with the addition of diurnal wind cycle at a rip-channeled beach. Modeled phytoplankton were released at two locations: offshore to examine how typical coastal phytoplankton taxa are transported onshore and increase  cell concentrations, and how surf zone taxa return to the surf zone; and inside the surf zone to examine how the coastal and surf zone phytoplankton taxa are retained in the surf zone and increase concentrations.

Study Site
Bathymetry and wave data were collected at Sand City beach (36 • 36 ′ 57 ′′ N, 121 • 51 ′ 15 ′′ W), Monterey Bay, California in the summer of 2010. This area receives a strong sea breeze during the daytime and little, if any, land breeze at night (Hendrickson and MacMahan, 2009), and we observed the same wind pattern during the fieldwork. Bathymetry was taken with a personal watercraft equipped with sonar and a Global Positioning System (GPS), and a person walking with a GPS for shallow bottom. This location is characterized as an intermediate beach (Wright and Short, 1984) and rip channels and shoals are developed (Figure 3). Wave data were collected by acoustic Doppler current profilers deployed at 11 m water depth. The average wave direction during the field experiment was perpendicular to the shoreline (Figure 3). Numerous studies described hydrodynamics of this beach (MacMahan et al., 2004Reniers et al., 2009Reniers et al., , 2010Fujimura et al., 2013Fujimura et al., , 2014Brown et al., 2015). Typically, bathymetric rip currents and corresponding eddies are formed here as described in Figure 2.

Hydrodynamic Model
Three-dimensional hydrodynamic model simulations have been performed with Delft 3D which comprises FLOW (Deltares, 2011a) and WAVE (Deltares, 2011b) modules. All hydrodynamic simulations followed Fujimura et al. (2014) except for modeling durations and wind and wave conditions. For all model cases, we used only regular normally incident waves based on the averaged measurement data, and included Stokes drift (Stokes, 1847), a wave-related time-averaged volumetric transport in the direction of wave propagation. Fujimura et al. (2014) showed that Stokes drift is an essential part of onshore particle transport.
We considered situations without wind and with a diurnal wind cycle. In the no-wind condition, the model was run only with wave-driven currents for 48 h. For the diurnal wind cycle, a simulation started with constant offshore wind (2 m s −1 ) for the first 12 h, linearly changed to onshore wind (6.0 m s −1 ) in 1 h and lasting for 24 h, then gradually returned to an offshore wind in 1 h. The wind again gradually shifted from offshore to onshore wind at 36 h, and the model was running with the constant onshore wind until 48 h. The first and third quarters of the simulation (0-12 h and 24-36 h) were considered as daytime, while the second and fourth quarters (12-24 h and 36-48 h) were simulated as nighttime.

Phytoplankton Transport Model
Modeled physical parameters (i.e., water currents, waves, eddy diffusivities and turbulence) and bathymetry were transferred to an individual-based model (IBM). We adapted the IBM code developed by Fujimura et al. (2014). We kept basic functions, such as the 4th order Runge-Kutta method for advection, random walk based on simulated diffusivities, and the boundary conditions from the open-source connectivity modeling system .
A basic vertical velocity of each phytoplankton cell was downward w p = −1.2 × 10 −5 m s −1 (1 m d −1 ), which is a typical sinking rate of phytoplankton (Smayda, 1970). As we mentioned earlier, sinking velocity of cells may be increased by turbulence; thus, we assumed that cells sink at w sink = −2.5 × 10 −3 m s −1 when the turbulence energy dissipation rate is ε > 10 −5 m 2 s −3 (Ruiz et al., 2004). Turbulence greater than this threshold can be seen in the surf zone and some part of bottom boundary layer.
The total velocity components u tot (cross-shore), v tot (alongshore), and w tot (vertical) are: where u adv , v adv , and w adv are advection flow velocities, and u diff , v diff , and w diff are random velocities of all three spatial components (cross-shore, alongshore, and vertical direction, respectively). u adv and v adv are Lagrangian velocities (u L , v L ), consist of Eulerian velocities (u E , v E ) and Stokes drift.
We also added a growth rate that increases under high turbulence condition (Savidge, 1981). The main purpose of the model is to examine the difference between phytoplankton concentrations in the surf zone and offshore, so the total number of cells in the system does not matter. Hence, in this model, not all phytoplankton cells grow, but only cells in high turbulence can grow. Every hour, randomly selected particles where ε > 10 −5 m 2 s −3 are doubled in cell number. This ε is the same value as the threshold for the turbulence-driven sinking. The number of randomly selected cells per hour is: where N t is total cells in turbulence, µ is a specific growth rate of phytoplankton and r is a fraction of increased growth rate due to turbulence. Here we use µ = 0.5 d −1 (Parsons et al., 1984) and r = 0.5 (Savidge, 1981). Separately, the DVM model proposed by Talbot et al. (1990) was applied to our model case with a diurnal wind cycle. In our model, DVM was considered only within the surf zone (X ≤ 100 m, see Figure 3). As suggested by Talbot et al. (1990), phytoplankton cells float on the surface during the day, and sink to the bottom at night. Floating cells are attached to air bubbles in the surf zone (Talbot and Bate, 1988b;Shanks et al., 2018). We assumed that the diameter of air bubbles was 100 µm, a typical size in the surf zone (Deane, 1997;Deane andStokes, 1999, 2002). Based on empirical data (Detsch, 1991), we estimated an upward velocity of the bubbles as w air = 6.0 × 10 −3 m s −1 . Sinking cells are assumed to be attached to sand grains in the water column (Talbot and Bate, 1988b). We chose a uniform grain size 0.3 mm in diameter, which is a typical size at Sand City beach (Gallagher et al., 2011) and can be suspended in the surf zone . A sinking velocity of the grain w sand = −1.45 × 10 −2 m s −1 was calculated with Zhiyao et al. (2008). Therefore, for the DVM model, Equation 3 becomes: w tot = w adv + w diff + w p + w air during the day w adv + w diff + w p + w sand during the night Frontiers in Marine Science | www.frontiersin.org In addition, we considered here phytoplankton growth during the day with vertical variation in the surf zone. For the DVM model, we assumed that cell division can be observed in almost all phytoplankton on the surface, half of the phytoplankton in the water column, and a very small portion of phytoplankton on the bottom (Talbot et al., 1990). When the DVM model included the phytoplankton growth in the surf zone, we applied a growth rate of 1.0 d −1 to cells at the depth of ≤ 0.25 m, 0.5 d −1 to cells in the water column, and 0 d −1 to cells at ≥ 0.25 m above the sea bed.
We tested the effects of increased sinking velocity and growth rate of phytoplankton owing to turbulence with nowind or diurnal wind condition, as well as the effect of DVM in the diurnal wind regime, on cross-shore phytoplankton transport. Each simulation run was 48 h. A total of 602 cells (86 alongshore × 7 vertical array) were released every hour from offshore (X = 410 m) to examine onshore transport of offshore phytoplankton taxa and a return rate of surf zone diatoms; or inside the surf zone (X = 60 m) to see an exit rate of offshore species and a retention rate of surf zone taxa. Only the last quarter of each simulation (36-48 h) was used for analysis. The earlier period (0-36 h) was used as a spin-up stage for initialization. The model cases and parameters are summarized in Table 1.

Cross-Shore Velocity Profiles
Alongshore-and time-averaged cross-shore velocity profiles are shown for the defined rip channels and shoals (Figure 3) within the surf zone (X < 100 m), and offshore (X > 200 m) in the total alongshore range (Figure 4). The result verified the theory that Stokes drift adds onshore forcing to Eulerian current. Shoal and rip current were apparent at X = 75 m, and bottom boundary streaming enhanced onshore bottom (X ≥ 200 m).
Onshore wind induced onshore surface currents in the daytime (Figure 4). At night, offshore wind altered the flow condition similar to that in the no-wind regime. Flow directions in the surf zone did not change with wind directions although magnitudes slightly changed. As expected, bottom streaming was suppressed during the onshore wind event because of the mass balance needed at the bottom to counter the onshore surface flow. Overall, the cross-shore current fields produced were consistent with the concept presented in Figure 1.

Transport and Distribution of Phytoplankton
Many offshore-released phytoplankton cells were ejected from the domain (Figure 5). Phytoplankton were not effectively carried in the surf zone without turbulence-driven sinking. Turbulence-driven growth had little effect as growth was only enhanced when cells were within the surf zone.
In the diurnal wind regime, the number of offshore phytoplankton that were transported toward shore was highest in the case that the cells had the turbulence-driven growth and sinking, followed by the case with DVM and depth-varying growth (Figure 5). DVM somewhat enhanced onshore transport, but it was more effective at increasing the concentration of phytoplankton within the surf zone when the growth function was added. Interestingly, some phytoplankton cells without Phytoplankton cells were released either at X = 410 m (offshore) or 60 m (within the surf zone). The wind condition was either no-wind or a diurnal wind cycle (daytime onshore breeze and nighttime offshore breeze). Diurnal vertical migration (DVM) was not included with turbulence-driven growth or sinking. Corresponding figures are indicated.
Frontiers in Marine Science | www.frontiersin.org 5 June 2018 | Volume 5 | Article 183 turbulence-driven sinking and DVM were also carried to shore ( Figure 5E). When phytoplankton were released in the surf zone, all cases showed cell concentrations higher in the surf zone than offshore (Figure 6). Sinking with sand grains enhanced retention of cells in the surf zone. Surf zone concentration of phytoplankton released within the surf zone was highest in the case of the cells with DVM and depth-varying growth ( Figure 6H); however, cell retention rate in the surf zone in this case was not high compared to the other cases as substantial quantities of phytoplankton seemed to be carried offshore ( Figure 6). Phytoplankton with turbulence-driven growth and sinking tended to be retained and subsequently increased in the surf zone ( Figure 6F). Similar to the offshore-released case, cells without any traits also tended to stay in the surf zone ( Figure 6E).
Cell concentrations in the surf zone would become much higher than offshore over time because onshore cell flux (Figure 5) was higher than offshore cell flux (Figure 6). All (except for Figures 5A,B) cases showed the common distribution trend of phytoplankton (per unit area), that is, highest cell concentrations were found in the rip channels, followed by Frontiers in Marine Science | www.frontiersin.org that over the shoals, and the smallest concentrations offshore (Figures 5, 6).

Comparison Between No-Wind and Diurnal Wind Case
Phytoplankton without sinking, growth, and DVM were dispersed vertically in the no-wind condition (Figure 7A), whereas those in the diurnal wind condition tended to be concentrated near the bottom and surface (Figures 7C,E).
Phytoplankton concentration at the surface extended offshore in the no-wind case (Figure 7A). A similar extension was formed at night in the diurnal wind case (Figure 7C), but not pronounced like in the no-wind case. Vertical velocities at X > 150 m during the day in the diurnal wind cycle directed downward (Figure 7F). Offshore vertical velocities (X > 200 m) at night tended to be slightly upward (Figure 7D), and neutral with a little variability in the no-wind regime ( Figure 7B).

DISCUSSION
This modeling study suggests that turbulence-driven sinking may be an important mechanism for transport of phytoplankton into the surf zone, just like the sinking behavior of competent larvae in turbulence (Fujimura et al., 2014). Once phytoplankton enter the surf zone, their growth may be enhanced by turbulence, resulting in high cell concentration in the surf zone. Our model showed that the diurnal wind cycle increased onshore transport (Figures 5A,D compared to Figures 5E,F) and reduced flushing rates of phytoplankton cells (Figures 6A,D  compared to Figures 6E,F). These results were probably due to onshore wind, the major difference between no-wind and diurnal wind regimes. The modeled phytoplankton cells had a relatively small fall velocity (w p = −1.2 × 10 −5 m s −1 ), which is 1-2 orders of magnitude smaller than the modeled vertical water velocities (Figure 7); thus, the vertical displacement of cells was affected by surrounding currents. In all cases, phytoplankton cells were brought to the surface by the upward flow in and near the surf zone (Figure 7). When wind was not applied, many of the cells were kept near the surface throughout the domain, and some cells were sinking or floating in the water column offshore (Figure 7A). In the diurnal wind cycle, cells floating near the surface were also carried offshore at night (Figure 7C), but upward currents kept them off from the water column (Figures 7C,D); furthermore, these cells were pushed back to the shore when the condition switched to onshore wind condition ( Figure 7E). During the day, downward currents in the water column enhanced downward velocities of phytoplankton (Figure 7F), and transported via the onshore bottom currents although these were weakened due to balancing with onshore surface currents (Figure 4). Therefore, wind-driven onshore currents were likely to enhance phytoplankton transport toward the shore, while the constant wave-driven flow regime without wind caused higher phytoplankton flushing rates. The results were consistent with a study by Hendrickson and MacMahan (2009) who indicated that diurnal sea breeze is important for cross-shore transport.
Other physical forcing can play a role as well. Stokes drift represented by difference between u L and u E (Figure 4) added onshore forcing. Without Stokes drift, onshore transport of materials at the surface may be slower as shown by Fujimura et al. (2014). Bottom onshore currents and benthic streaming owing to wave forcing are also important especially when plankton are in the wave boundary layer. Alongshore variability (rip channels and shoals) can be another factor that increases onshore phytoplankton transport. Bathymetrically-controlled rip currents enhance shoaling to balance cross-shore currents. This bathymetric feature often forms a rip circulation with an eddy at its center where plankton may be trapped (Talbot and Bate, 1987;Fujimura et al., 2014;Shanks et al., 2017). Furthermore, flatter beach slopes (i.e., more dissipative beach) should be more conducive to onshore transport of plankton (Shanks et al., 2010(Shanks et al., , 2018Fujimura et al., 2013;Morgan et al., 2016Morgan et al., , 2017. Our data showed that the concentrations of phytoplankton were highest in the rip channels, followed by that over the Frontiers in Marine Science | www.frontiersin.org shoals, and lowest offshore. This pattern is the same as the larval distribution at the same beach (Fujimura et al., 2014;Morgan et al., 2016Morgan et al., , 2017. High concentrations of plankton in the rip channels are possibly a result of the convergence of feeder currents (Figure 2). Our results partially agreed with Shanks et al. (2018) that phytoplankton concentrations in the surf zone were much higher than offshore. However, contrary to our study, their data showed that cell concentrations over the shoals were lower than offshore. This discrepancy might be due to the sampling scheme. They collected samples at 1 m depth in the surf zone and throughout the water column offshore, while we averaged phytoplankton cells in the water column everywhere. Hence, we recalculated the time-averaged number of cells over the shoals at 1 m ± 0.5 m. Nonetheless, modeled phytoplankton concentrations were still higher over the shoals than offshore in all cases although shoal to offshore ratios became smaller (e.g., for DVM with depth-varying growth case, 12:1 with original calculation, 2:1 with recalculation). More plausible reason for the discrepancy between Shanks et al. (2018) and our result was a lack of significant environmental parameters and/or phytoplankton traits and further investigation will be necessary for the model to reproduce such a distribution pattern of phytoplankton.

Limitations and Generality
Phytoplankton cells are suspended occasionally (e.g., in the ocean surface mixed layer) which has been supported by some studies showing that downward velocities of phytoplankton slowed in turbulent waters (e.g., Ruiz et al., 1996;Deleersnijder et al., 2006); nevertheless, this phenomenon can be seen only when turbulence level varies vertically (Ross, 2006;Macías et al., 2013). As the surf zone in our model was a more or less homogeneous layer with the turbulence energy dissipation rate ε > 10 −5 m 2 s −3 , the value used as the threshold obtained from the empirical data (Ruiz et al., 2004), the reduction of vertical velocity was ignored here. Vertical motion in turbulence is species-specific (Ruiz et al., 2004), while most studies on vertical velocities in turbulence were not conducted with real phytoplankton cells, so more observational and laboratory data will be valuable for more accurate models.
Turbulence may enhance cell growth by increasing nutrient uptake and/or exposure to light. Sullivan et al. (2003) reported that growth rates of dinoflagellate species increased with turbulence levels up to ε ≈ 10 −4 m 2 s −3 , but decreased at ε ≈ 10 −3 m 2 s −3 although they still had fairly high growth rates. If we consider their result, the model may be improved by including a growth function that varies with cross-shore varying turbulence energy dissipation rates in the range of 10 0 -10 −5 m 2 s −3 (shoreline-surf zone edge). Also, the effect of turbulence on growth rate is species-specific as turbulence influences the growth rate of some species negatively (Peters and Marrasé, 2000), while others grow faster in turbulence (Davis et al., 1953;Savidge, 1981;Hondzo and Wüest, 2009). Moreover, the surf zone often receives nutrient-rich freshwater from land through the sand of the beach, and there may be spatial nutrient variability (i.e., the closer to shore, the higher in nutrients). In addition, the surf zone is a relatively shallow region, so generally phytoplankton are exposed to high irradiance, resulting in fast growth rates.
Simulations of phytoplankton concentration can be improved by developing cell growth model with corresponding environmental variability.
Our DVM with depth-varying growth model was based on Talbot et al. (1990). This model case showed the most successful onshore phytoplankton migration in our simulations; however, it is not known whether this mechanism applies to any species. As mentioned previously, the ability to float by attaching to air bubbles has been observed in surf zone diatom species (Lewin and Schaefer, 1983), but offshore phytoplankton can also attach to bubbles. In fact, Shanks et al. (2018) collected foam from the sea surface in the rip current, and found that phytoplankton are often highly concentrated in the foam. Those phytoplankton consist of a lot of offshore diatoms and dinoflagellates, and a tiny percentage of surf zone species. It would have been interesting if they had determined the vertical distribution of phytoplankton and sampled at night like Talbot and Bate (1988a) performed.
The model can be improved in many aspects. It is necessary to collect species-specific information including habitat, growth rate, sinking and floating rates, cell size, and morphology. These parameters often interact with each other and the modeling is not trivial even for single species. If we consider that sand grains are sinking agent of phytoplankton, their size variability in the surf zone from observed data (e.g., Gallagher et al., 2011) and their dynamics from numerical model (e.g., Reniers et al., 2013) would help to improve the phytoplankton biophysical Lagrangian model. Likewise, if the target species interact with air bubbles, dynamics and size distributions of air bubbles produced by breaking waves in the surf zone from empirical data (e.g., Deane, 1997;Deane andStokes, 1999, 2002) or numerical simulations (e.g., Ma et al., 2011) could be incorporated into the biophysical model. Air bubbles may be entrained by "obliquely descending eddies" (Nadaoka et al., 1989), which our model did not include. Another missing surf zone process in our model is a breaking wave roller (Feddersen, 2007;Reniers et al., 2013) that may entrain phytoplankton cells and transport them toward the shore. Moreover, wave conditions also influence onshore transport of plankton. Concentrations of phytoplankton in the surf zone at Sand City beach were higher between 0.5 and 1.0 m wave heights (Shanks et al., 2018). Most zooplankton concentrations in the surf zone tended to be negatively correlated with wave heights at Sand City beach  and at a more reflective beach (Shanks et al., 2015;Fujimura et al., 2017;Morgan et al., 2017). Wave spectra may need to be applied for more realistic situations. Random wave groups generate infragravity waves and surf zone eddies that may add extra forcing to onshore transport of plankton and material within the surf zone (MacMahan et al., 2004;Reniers et al., 2010;Fujimura et al., 2014).

CONCLUSIONS
We found that both biological and physical parameters controlled onshore delivery of phytoplankton and the concentration of the cells in the surf zone. Both wave and wind forcing play an important role in phytoplankton transport. More modeled offshore phytoplankton cells tended to reach the surf zone in the diurnal wind cycle regime than that in the no-wind condition. The cells were highly concentrated in the rip channels by current convergence and rip circulations. Turbulence-induced growth and sinking enhanced the number of phytoplankton cells in the surf zone. DVM also facilitated onshore phytoplankton transport, and depth-varying growth increased the cell concentrations in the surf zone.
The model parameters could be improved with additional empirical data. Further investigation on dynamics of waves, currents, air bubbles, and sand grain will help to develop the model. Species-specific morphological, physiological, and behavioral characteristics of phytoplankton need to be studied as well. Once the model with necessary parameters is validated, it may be coupled with other ecological model such as NPZD and nested in a larger scale for more practical applications, including HAB predictions and fisheries management.

AUTHOR CONTRIBUTIONS
Designed the investigations: AF, AR, AS, JM, SM. Developed the model: AF, AR, CP. Performed the simulations and data analysis: AF. Prepared manuscript: All authors.