Numerical Modeling of Microbial Fate and Transport in Natural Waters: Review and Implications for Normal and Extreme Storm Events

: Degradation of water quality in recreational areas can be a substantial public health concern. Modelscanhelpbeachmanagersmakecontemporaneousdecisionstoprotectpublichealthatrecreational areas, via the use of microbial fate and transport simulation. Approaches to modeling microbial fate and transport vary widely in response to local hydrometeorological contexts, but many parameterizations include terms for base mortality, solar inactivation, and sedimentation of microbial contaminants. Models using these parameterizations can predict up to 87% of variation in observed microbial concentrations in nearshore water, with root mean squared errors ranging from 0.41 to 5.37 log 10 Colony Forming Units (CFU) 100 mL − 1 . This indicates that some models predict microbial fate and transport more reliably than others and that there remains room for model improvement across the board. Model refinement will be integral to microbial fate and transport simulation in the face of less readily observable processes affecting water quality in nearshore areas. Management of contamination phenomena such as the release of storm-associated river plumes and the exchange of contaminants between water and sand at the beach can benefit greatly from optimized fate and transport modeling in the absence of directly observable data.


Introduction
Water systems and the recreational opportunities that they afford bring millions of people outside each year, especially during warm weather. Over 75% of people traveling in the summer visit beaches, and in Chicago alone 20 million people go to Lake Michigan beaches annually, on average [1,2]. To protect public health and the safety of beachgoers, the Beaches Environmental Assessment and Coastal Health (BEACH) Act of 2000 requires routine monitoring of coastal water quality at both marine and freshwater beaches across the USA [3]. This monitoring, however, often involves obtaining samples and either culturing for fecal indicator organisms (FIO) such as E. coli or enterococci or using quantitative Polymerase Chain Reactions (qPCR) to determine FIO concentrations in the water. These approaches take time, leading to a delay of up to 24 hours before obtaining water quality information to effectively manage beach usage for public health. This delay can be the difference between keeping beachgoers safe by advising against beach access and putting them in danger by keeping a contaminated beach open for recreation.
Numerical mechanistic models have shown variable success in predicting microbial concentrations in coastal and beach systems across space and time. As knowledge of aquatic systems and their various influences on waterborne microorganisms has progressed in recent years, predictive capacity of models has increased as well [13,16,17,21]. However, with a still incomplete knowledge of how components of the aquatic environment interact to influence microbial fate and transport, improving the representation of processes as well as of source behavior, parameter identification and model evaluation remains an evolving process.
This process is further complicated by the often rapidly changing conditions within aquatic environments. For example, Lakes Michigan and Huron in the Laurentian Great Lakes, USA, have become significantly clearer since the 1990's, in response to an invasion by dreissenid mussels [22][23][24]. This clarification likely impacts microbial survival in the water, due to the subsequent changes in sunlight extinction and solar microbial inactivation rates in the water [25]. Similarly, changes to microbial sources and hydrodynamics associated with climate change [26][27][28][29] may lead to changes in how microbial fate and transport can be effectively modeled. Not only is sea level predicted to rise due to climate change, but frequency and intensity of storm events are projected to increase for many regions as well [30]. Sea level rise will undoubtedly change the layout of beach areas, leading to phenomena such as shoreline erosion and increased foreshore areas that are susceptible to microbial transfer between water and sand [31]. At the same time, a predicted increase in the frequency and intensity of storm events will likely lead to increased runoff of urban and agricultural contaminants to rivers, which will in turn flow to coastal areas and potentially impact beach water quality [32].
Herein, we aim to compare the various approaches to modeling microbial (i.e., FIO) fate and transport in coastal aquatic environments using coupled mechanistic hydrodynamic and transport modeling approaches, to determine how these approaches impact model predictive capability. Examples of such models include the Finite Volume Community Ocean Model (FVCOM, [33]), the Princeton Ocean Model (POM, [15,20,34,35]), the Aquatic Ecosystem Model 3D (AEM3D, [17]), Delft3D [36], and the Environmental Fluid Dynamics Code (EFDC, [37]) to name only a few. For the purposes of this review, we will concentrate on these mechanistic and process-based models of FIO fate and transport, but statistical and data-based modeling approaches are briefly discussed to the extent that they can support mechanistic modeling efforts. We then discuss the application of such approaches to emerging modeling questions surrounding the simulation of storm-associated river plumes and microbial exchange between beach sand and water. While nearshore environments harbor a variety of microorganisms, including bacteria, viruses and fungi, this work will focus on modeling of FIO such as E. coli, enterococci, and coliforms. Understanding how numerical models predict water quality will lend insight into which approaches may be most appropriate for modeling recreational water quality to ensure public health in the face of climatic and environmental changes.

Modeling Hydrodynamics and Microbial Fate & Transport
Numerical modeling of hydrodynamics in coastal environments depends heavily on the physics of the water body, meteorological conditions over time and the impacts of river/estuary inputs. Water moves in three spatial dimensions and over time, so utilizing a three-dimensional hydrodynamic model is key to adequately simulating water movementin coastal areas. Unstructured grid models (e.g., FVCOM) may have advantages in simulating coastal water quality due to their ability to accurately represent nearshore features such as irregular coastlines, barrier islands and sandbars, harbors, breakwaters etc.
In these equations, x, y, and z represent the east, north and vertical directions, while u, v, and w denote the velocity components in the x, y, and z directions, respectively (m s −1 ). fu, fv are the Coriolis terms, F u , F v are the horizontal diffusion terms in the x and y directions, respectively and F w is a vertical diffusion term (m 2 s −1 ). K m denotes a vertical eddy viscosity coefficient (m 2 s −1 ) and ρ 0 represents the density of water (kg m −3 ). Air pressure at the water surface is denoted by p a , hydrostatic pressure is represented by p H , and q is the non-hydrostatic pressure (all in Pa) [33]. Using these equations, models can effectively account for the effects of temperature, density, and Coriolis force on water movement over time. The effects of waves in the nearshore environment (e.g., wave-current interactions, bottom shear stress) can be simulated using coupled hydrodynamic and spectral wave models such as the FVCOM-Surface Wave (FVCOM-SWAVE) model [40].
Building on the general hydrodynamic model, microbial fate and transport associated with diffusion, dispersion, advection, and mortality within an aquatic system can be simulated [9,12,13,[15][16][17]20,21,41]. The governing equation for microbial fate and transport is based on the advection-dispersion-reaction (ADR) equation, formulated in terms of FIO concentration (Equation (5)). This equation includes terms for advection, diffusion/dispersion in the water column, and microbial decay [16,33].
In the equation, C corresponds to the microbial concentration (Colony Forming Units (CFU) 100 mL −1 ), k represents the overall microbial decay rate (d −1 ), and K V and K H are the vertical and horizontal mixing coefficients, respectively (m 2 s −1 ). Horizontal and vertical mixing are described using the Smagorinsky and Mellor-Yamada 2.5 level turbulence parameterizations [33,39]. Microbial decay can depend on factors such as microbial taxon and base mortality, water temperature and chemistry, attachment to and detachment from suspended solids, settling after attachment to suspended solids, sunlight inactivation, and interactions with other biota in the aquatic environment [21]. Because of its dependence on these factors, the decay term in mechanistic models has taken on many forms within the literature and is often a combination of multiple terms.

Boundary and Initial Conditions
Key drivers of hydrodynamics in the nearshore region include winds and/or tides [42,43] and riverine/estuarine flows [9,[44][45][46]. These are often governed by local conditions such as bathymetry, wind stress, Coriolis force, and water temperature, which can vary spatiotemporally. For example, seasonality of thermal or density stratification in large lakes such as Lake Michigan can impact circulation as well as buoyancy of contaminant plumes, leading to differential impacts on water quality and hydrodynamics throughout the year [47,48]. Similarly, estuarine exchange flows and the corresponding changes in vertical density stratification (and hence vertical mixing) are controlled by the along-channel wind component which changes seasonally [49]. In addition to seasonal-temporal variability in fluid properties, there are spatial influences that control hydrodynamics. The impacts of Coriolis force can vary by latitude as well as water body size, with effects becoming significant for large (>5 km width) lakes and at high latitudes [50]. Because these effects can vary spatiotemporally while also influencing large-scale hydrodynamics by influencing stratification, it is important to specify related model boundary conditions as realistically as possible.
Boundary conditions for the momentum equations are well-known and include wind stress on the surface of the water column and bed friction on the lake/seabed; additional details are available in Chen et al. [33,39]. For the FIO transport model, the nature of the source(s) dictates the type of boundary conditions used. For beaches impacted by riverine sources, monitoring data collected at the river mouth can be used to provide boundary forcing for the FIO transport model. However, most mechanistic models use small time steps (on the order of seconds to minutes) while monitoring data are collected less frequently (e.g., weekly or bi-weekly), introducing significant uncertainty into the modeling due to a mismatch between the model time steps and forcing data. In addition, model inputs are generally specified at regular intervals while monitoring data may be irregular and with gaps (e.g., missing data during weekends).
One way to address this limitation is to use calibrated watershed models [51][52][53] or statistical models to generate high-resolution boundary forcing data at the river mouth when high-resolution discharge data are available (e.g., at a United States Geological Survey (USGS) gauging station in the USA). Several researchers have exploited a potential correlation between river discharge (Q) and FIO (e.g., E. coli) concentrations (C) and have used statistical relations between Q and C to generate high-resolution tributary loading data for nearshore FIO models [16,17,20]. Compared to the use of sparse monitoring data to represent tributary loading in FIO models, these approaches have promise, as they can better describe rapid changes in loading and may be suitable for simulating the impacts of extreme storm events on microbial water quality in coastal areas. For beaches with no known riverine sources, observed FIO dynamics may be driven by local sources including birds [12,54,55], resuspension of bottom sediment-bound FIO [45,46] and shoreline sand [7,56]. A comprehensive analysis of the relative importance of the different sources calls for detailed modeling of hydrodynamics including currents and waves, sediment transport, sediment-FIO interactions and relatively fine computational grids to capture the impact of shoreline birds and sand. Because initial conditions for hydrodynamic models of lakes and reservoirs may specify a waterbody that is initially at rest (zero velocity components), a spin-up period is often used to allow models to catch up with observed data. For FIO models, the initial concentration of FIO may include a small non-zero background value and model spin-up time may allow for increased water quality predictability over the model simulation period.

Components of the Microbial Decay Function
One form of the microbial decay function from Equation (5) is represented by Equation (6), where total microbial loss (k) is characterized by the combination of a base mortality rate (k b1 ), light inactivation rate (k bi ) and settling rate (k bs ) [41,57]. All of the decay terms in the equation are daily decay rates (units of d −1 ). k = k b1 + k bi + k bs (6) The factors that influence these three terms, though, can be variable and specific to the model domain, context, and aims.
Liu et al. [13] further separated this basic decay function for Lake Michigan, yielding a function that accounted for bacterial loss due to settling and light inactivation, with a temperature correction factor to justify changes due to temperature variations from 20 • C (Equation (7)). In this model, f p is the fraction of FIO particles attached to suspended solids (unitless), v s represents settling velocity (m d −1 ), H is the depth of the water column (m), k I is the inactivation rate associated with solar radiation (m 2 W −1 d −1 ), I t denotes the solar irradiance at the surface of the water at time t (W m −2 ), θ is the temperature correction factor (unitless), and T is the water temperature ( • C).
The right-hand side of this equation is composed of a sedimentation/settling term, a light inactivation term, and a temperature correction factor, when read from left to right. Though it incorporates terms for temperature, settling, and sunlight impacts on FIO, this form of the decay function fails to account for other potentially important influences.
Microbial survival in aquatic systems can be subject to impacts due to temperature, salinity, light penetration and inactivation, predation, competition for resources, nutrient availability, and other natural (or base) mortality. Similarly, conditions such as the presence of aquatic vegetation and dreissenid mussels in a water body can impact bottom shear stress, thereby influencing sedimentation and resuspension of FIO. These factors can be difficult to characterize, especially when they vary between systems. In marine or brackish ecosystems where salinity may vary over time and space, it can have a substantial impact on the survival and persistence of several microorganisms [21,[58][59][60][61][62][63][64][65][66]. Similarly, in systems with high levels of mixing or turbidity, FIO may not settle out of the water column as much as they would for less well-mixed systems that would foster settling [67]. Finally, light inactivation has been shown to play a significant role in FIO decay via inactivation in natural waters [68,69]. This is especially true for oligotrophic or low-turbidity waters that do not have high levels of suspended solids that can serve as refugia for FIO [25]. Impacts can also vary between different microorganisms. Cabelli [70], Colford et al. [71], and Schang et al. [72] document differential sensitivities to environmental pressures between E. coli, Cryptosporidium parvum, Giardia lamblia, and Campylobacter spp. in natural waters. Sanders et al. [73] tested the impact of organism-dependent sensitivities to environmental survival influences on model predictive ability. They found that a model with the same environmental influences underpredicted E. coli concentrations, but overpredicted concentrations of both total coliforms and enterococci, compared to observations.

Dark Mortality, Base Mortality and Temperature/Salinity Dependence
Base mortality (k b1 ) and dark mortality (k d ) terms are often used in models to account for the natural decay that microbes undergo, independent of sunlight effects. The value of k d can be highly variable, depending on the geographic location and the microorganism of focus [15,16,21,[74][75][76][77]. Models have been presented using k d values as low as 8.6 × 10 −5 d −1 [15] and in vitro experiments have yielded k d rates of up to 2.2 d −1 [74]. Once k d is established, it is then possible to correct for other factors affecting k b1 , such as temperature, salinity, or pH of the water.
Dark mortality rates are developed for a reference water temperature of 20 • C, so an adjustment to account for variability in temperature is also needed. Microorganism base mortality terms are adjusted for temperature using the Arrhenius equation [15,18,78]. Resulting temperature correction factor values can range from 1.04 to 1.11 [78], but frequently are assumed to be 1.07 [79]. This adjustment indicates a strong temperature dependence, with a doubling of the mortality rate for every 10 • C increase in temperature [41]. The resulting full formulation of the base mortality term is thus represented by equations 10 and 11, where θ represents the temperature correction factor.
A majority of existing models use some version of this formulation to determine base mortality rate of microorganisms in natural waters (Table 1). Notable exceptions were found in models from Liu et al. [80], McCorquodale et al. [81], Hipsey et al. [21], Rehmann and Soupir [78], Servais et al. [82,83], and de Brauwere et al. [84]. Rather than using a dark mortality rate to calculate base mortality, Liu et al. [80] use the time to inactivate 90% of microorganisms in the dark (t 90 ). McCorquodale et al. [81] use a curve-fitting procedure on field-collected data to determine the impact of salinity on base mortality. Hipsey et al. [21] and Rehmann and Soupir [78] incorporate the effects of salinity and pH on mortality (c S M and c pH M , respectively), sensitivity of the microorganism to salinity and pH (β and K δ pH M , respectively), nutrient limitation (f LIM ) and dissolved organic carbon concentration (DOC L ) when calculating base mortality. Following Servais et al. [82,83], de Brauwere et al. [84] use a logistic relationship between microbial decay and temperature to determine k b1 .

Solar Inactivation Terms
It is generally accepted that incoming solar radiation affects the survival of microbes in water systems [69,85,92,93]. This is especially true in clear, oligotrophic waters, where solar inactivation can be a predominant influence on microbial survival [25,68]. Many mechanistic fate and transport modelers recognize the impacts of solar inactivation on microbial survival in water and include inactivation parameters in their models.
Accounting for solar irradiation in natural waters inherently involves the calculation of the light extinction rate within the water column. The amount of light penetrating the water column declines exponentially with depth and is influenced by the turbidity or clarity of the water, such that clearer water yields a lower light extinction rate than turbid water. Lower light extinction rates, in turn, yield more intense solar radiation at deeper depths in the water, leading to higher microbial solar inactivation rates [41].
Nearly all of the models which account for solar inactivation use some variation of Equation (7) in which the Beer-Lambert Law (Equation (12), [94]) is used to model the variation of solar radiation with depth. In the Beer-Lambert equation, z is vertical coordinate of depth (m) and I z represents the amount of solar radiation at vertical coordinate z (W m −2 ). I 0 is solar radiation at the water surface (W m −2 ), and k e is the light extinction rate (m −1 ) [95]. It is important to distinguish between k I from Equation (7) and k e in Equation (12). In Equation (7), the k I term is an inactivation rate for FIO as a result of solar radiation [16,73,80], whereas k e in Equation (12) is the rate of light extinction with depth in the water column [21,41,95].
Model type can have a significant impact on the variables used in parameterization of solar inactivation effects on FIO. Models may employ either total depth (H) or the vertical coordinate (z) within their solar radiation parameterizations, depending on the model context. In two-dimensional model frameworks, conditions within the water column are often vertically-integrated. For these cases, fate and transport models use a single depth variable (H) and a single solar radiation variable (I t ) to account for potential variability in the vertical dimension [13,41,45,78]. Fully three-dimensional models, in contrast, explicitly define conditions at different depths in the water column via their incorporation of the vertical coordinate variable z within their parameters. For example, fully three-dimensional models will often incorporate variables for solar radiation at the water surface (I 0 ) and solar radiation at depth z (I z ) to capture differences with depth in the water column [15,16,21,35,59,74,85,86,90,91,96].
A variety of approaches have been used for characterizing the effects of solar inactivation on FIO fate and transport. In some cases, models do not include solar inactivation terms at all [17,51,76,77,84,88], often because the water is so turbid that solar effects are assumed negligible compared to other environmental influences. Others use either a microbial decay rate solely as a function of incoming solar radiation (Equation (7)) or as a function of the light extinction rate and depth in the water (Equation (12), Table 2). Hipsey et al. [21] expanded the description of solar inactivation effects on FIO in their generic modeling framework, including specific terms for dissolved oxygen (DO), pH, salinity (S), and solar bandwidth (b). Garcia-Alba et al. [18] included terms corresponding to day length (DL) and fraction of solar irradiance that is in the UV spectrum (f UV ) as well as the typical light extinction rate, solar inactivation rate, solar irradiance and depth terms seen in other models.

Sedimentation Terms
In addition to solar inactivation, attachment to suspended solids and settling out of the water column is another significant driver of FIO losses in aquatic environments. 80-100% of total coliforms and E. coli have been shown to readily attach to suspended particles in the water column [100], and viruses have also been shown to easily attach to particulate matter and settle out of suspension [101].
Similar to the solar inactivation term, several published models do not incorporate sedimentation effects on microbial fate and transport [44,73,87,88,90,98,99,102,103]. In models that do incorporate sedimentation losses, settling terms most frequently use parameters representing settling velocity (v s , as calculated using Stokes' Law), vertical coordinate (z) or the total water column depth (H), and the fraction of the FIO that is attached to particles (f p , Table 3). In many cases, sedimentation terms are also subject to temperature correction, in the same manner that base mortality and solar inactivation terms utilize temperature correction factors [13,16,76,77], to acknowledge the fact that overall loss of FIO increases with temperature.  In their generalized sedimentation term, Hipsey et al. [21] expanded upon the simplified sedimentation terms used in most other models. This expansion accounts for various particle size classes (N s ), particle and attachment surface areas (A p and A s , respectively), and settling velocities for attached (v s ) and unattached (v c ) FIO.
Bravo et al. [20] and Thupaki et al. [35] incorporated sedimentation effects by including them in the vertical advection term of the 3D ADR equation (Equation (5)). As a result, the ADR presented is Equation (13) and the microbial decay function (kC) only includes terms for base mortality and solar inactivation in these models.

Model Testing and Evaluation
There is a large number of processes influencing FIO fate and transport and it can be difficult to identify a correct conceptual model that acknowledges process interdependencies over wide ranges of environmental variables of interest. Therefore, it is difficult to fully test FIO models across environments, and within the same environment, across different time periods (e.g., dry vs. wet weather events, "normal" vs. extreme events). While calibrated FIO fate and transport models have the potential to aid management by providing near real-time predictions, a majority of the published papers report results of model back-testing (or history matching, see Bredehoeft and Konikow [104]).
To evaluate the goodness of fit between models and observational data as well as to identify superior model formulations (by comparing different models), the use of multiple model evaluation metrics may be more beneficial [104,105] than the use of a single metric such as the coefficient of determination (R 2 ) or root mean squared error (RMSE). This is due to the fact that no single model performance metric captures all aspects of the data and simulation results, and all metrics have known limitations. In the context of FIO and beach management, evaluating models using the confusion matrix and concepts of sensitivity and specificity [19,106,107] have proven to be useful, especially for the practical application of issuing beach advisories and closings.
Existing, published models have been tested in a number of ways. Most model testing protocols, particularly those in more recent modeling studies, involve statistical analysis of comparability of model results to observed data. A majority of published models have used RMSE or R 2 as model performance metrics (Table 4). Other statistics such as Normalized RMSE (NRMSE), Mean Absolute Error (MAE), Nash-Sutcliffe Efficiency (NSE) [108], Percent Bias (PBIAS) [109] or the Refined Willmott Index of Agreement [110] have also been used by Boye et al. [90], Liu et al. [77] and Feng et al. [98]. A subset of the published models was qualitatively assessed, often with comparisons to other models replacing the more quantitative RMSE, R 2 , MAE, NSE and Wilmott index statistics [21,38,59,76,78,84,85,87,88]. In all applicable cases, the RMSE or MAE values have units of log 10 FIO CFU 100 mL −1 , while NRMSE units are percentages and NSE values are unitless. Good agreement with coliform mortality rates.
Mancini [59] Freshwater Lake/River Coliform Qualitative Evaluation Model was based on empirical relationships from lab/field data.
Auer and Niehaus [85] Freshwater Lake/River Fecal Coliform Qualitative Evaluation Model output is comparable and consistent with observed bacterial loads during wet and dry weather events.

Qualitative Evaluation
Fairly good prediction of observed microbial concentrations, but a general underestimation by the model.  Based on R 2 to evaluate model performance, the model of Garcia-Alba et al. [18] produced one of the best descriptions of observed data among the models considered here (R 2 = 0.87). Their model incorporated temperature-and salinity-dependent base mortality (k b1 = k d + k salinity θ T−20 ) and solar inactivation terms accounting for day length and fraction of irradiance composed of UV radiation ). Based on RMSE alone, one of the published models that best approximates observed microbial concentrations is detailed in Thupaki et al. [15] (RMSE = 0.41 log 10 CFU 100 mL −1 ). Within this model, k b1 = k d θ T−20 , k bi = (k I I t )·θ T−20 , and k bs = ∂( f p v s C) ∂z θ T−20 . Although comparison of these model frameworks can lend insight into which ones may best simulate microbial water quality, it is important to note that the models were developed under varying contexts. One published approach attempted to develop a generic water quality model, to be used across environments and target microorganisms [21]. This model framework led to complex terms within the decay function, often including parameterization for salinity, pH, dissolved organic carbon concentration, varying particle sizes and settling velocities, and variable sensitivity of the microorganism to such environmental changes. The resulting validation of the model indicated that the generic model did not outperform other existing models in its prediction of contaminant fate and transport. Despite this lack of substantial improvement over existing models, the generality of this model may be attractive to researchers looking for a single model to predict water quality under various conditions.
A lack of ancillary data may provide a confounding factor in the use of generic models such as the one described in Hipsey et al. [21]. In many cases, models are developed without the use of DOC, pH, temperature and salinity data, instead relying on hydrodynamics and meteorology to model FIO fate and transport. Likewise, additional water quality data such as DOC, pH, temperature and salinity are often not collected or available for model development, potentially hindering the usage and applicability of such a generic model. This could, however, indicate not that generic models may be less useful than specific and localized models, but rather that in situ DOC, pH, temperature and salinity data should be collected as part of the water quality monitoring process. For example, while the focus of most FIO modeling efforts is to reproduce the observed FIO concentrations, if no temperature data are collected in the nearshore region the ability of the coupled FIO-temperature-hydrodynamic model to accurately represent FIO decay is questionable as base mortality, solar inactivation, and sedimentation are often functions of temperature. A majority of the models reported in the literature use microbial decay formulations with a solar inactivation term and use FIO data collected during the daytime. Previous research shows that the highest levels of FIO are typically observed during the early morning hours (e.g., 6:00 AM) [111] due to the absence of solar radiation the previous night. Therefore, modeling the nighttime variation of FIO is important to correctly describe FIO levels in the morning [112]; however, since monitoring data are not collected at night, this aspect has not received much attention in the modeling literature.
In the absence of the specific data needed for the generic water quality model described above, selection of an appropriate modeling framework should be based on the target FIO as well as the environmental and hydrological context of the model.

Applying Microbial Fate and Transport Models to Extreme Storm Events
Numerical simulation of microbial water quality has evolved in recent years, as the dynamics of processes such as solar inactivation have become clearer. Even so, the incomplete knowledge of influences within aquatic systems on water quality indicates that there is still room for model improvement. Likewise, climate and land use/urbanization changes provide additional contexts for the prediction of microbial water quality [26]. Although the link between extreme precipitation and waterborne disease outbreaks is well known [29,113], the current generation of FIO models can be further refined and tested for their ability to reproduce observed dynamics during extreme storm events. Major areas of research impacting coastal water quality from the perspective of extreme storm events may include the exchange of FIOs between water and sand at beaches, the fate and transport of FIOs in storm-associated river plumes and the expansion of water quality monitoring research into microbial source tracking and environmental DNA (eDNA) for use in public health contexts.
The interaction between water and sand at the beach, and its impacts on recreational safety and water quality, has been an active area of discussion in recent years [56,[114][115][116][117][118][119][120]. Microorganisms in beach sands have been cited as potential sources of contamination and swimmer infection as early as 2003 [115]. The microbial community within beach sands is unique [15] in that it can serve as either a sink or a source of FIO to the adjacent recreational water, depending upon wave energies, currents and the movement of the water. When wave energy is low, FIO often get deposited from the water and into shoreline sand where they can form biofilm communities, while higher wave energy frequently leads to the release and re-suspension of FIO from the shoreline sand into the water [56,116]. These sand-based sources and sinks can heavily impact spatial and temporal trends in FIO concentrations at beaches. Further, the Intergovernmental Panel on Climate Change (IPCC) has predicted increases in wind speeds and wave heights/energies in mid-and upper latitudes as a result of climate change [121]. Similarly, the IPCC has predicted sea level rise in coming decades, a phenomenon already being observed, leading to changes in the beach face and the intertidal zone that is impacted by wave deposition/resuspension of FIO [56,122]. This will likely lead to increases in wave-induced FIO release from sands and into recreational water. Because of the potential for climate change to significantly impact sand-water exchange of FIO at beaches, it will be integral for numerical models to include sand-sediment-water interactions when predicting microbial water quality. Currently there are gaps in our understanding of these sand-sediment-FIO related processes and there is a need to further refine our mechanistic modeling approaches based on high-quality field observations and datasets which are often lacking. This will be especially important in substantially wave-impacted beach areas, to improve upon model predictions that exclude sand/sediment parameters [35,46].
Such models accounting for sand-sediment-water interactions at the beach may take inspiration from modeling frameworks that incorporate sedimentation. For instance, the modeling approach developed by Hipsey et al. [21] includes terms for various particle size classes, accounting for differential resuspension effects on "fine" and "coarse" particles. Fine particles require lower bed shear stress values for resuspension, compared to coarse particles, so it may be important to differentiate between the readily resuspended particles and those that are less likely to resuspend after deposition [98,100,123].
After simulating sediment transport as a function of particle size, sediment-FIO interactions can be modeled using attachment-detachment kinetics following those established for subsurface transport models [124].
An additional concern related to how climate change will impact recreational water quality involves storm-and runoff-associated FIO at coastal areas. For many regions, including mid-latitudinal coastal areas, climate change is expected to lead to increasingly frequent and intense storms [121]. Not only will these intense storms make recreation at beaches dangerous via rip tides, rip currents and strong waves, but they will also send increased volumes of potentially contaminated runoff and river water downstream, to be released to coastal areas [125]. As a result, recreational beaches may be expected to experience the impacts of more frequent and larger storm-associated river FIO plumes.
Effective prediction of the coastal water quality impacts from river FIO plumes will be helpful in not only understanding an additional source of contamination to recreational areas but will also aid in the management of beach resources for public and environmental health. This simulation will require the extension of existing numerical modeling approaches to include the determination of FIO concentrations in dynamic river plumes as well as reliable estimation of plume dynamics.
A number of studies have reported a "first-flush" effect for FIO (e.g., [122,123]) in which elevated FIO levels were observed following storm events with levels declining in later portions of the storm event and in subsequent events over a season. However, other researchers did not report such an effect [126]. These differences can be attributed to different runoff characteristics of watershed areas, so linking coastal water quality models with well-tested watershed models of FIO is expected to help address current limitations of nearshore FIO models [127]. For example, microbial composition and concentrations in runoff depend on upstream land uses; runoff from rural/agricultural watersheds is likely to have different water quality concerns than runoff from urbanized or forested catchments [128][129][130]. These differences are magnified during first flush phenomena and heavy storm events, where FIO can be released from soils and into the water, leading to high FIO loads in rivers that can then degrade coastal water quality. Therefore, calibrated upstream watershed models can be beneficial for modeling of coastal water quality in response to storm-associated river plume releases, simply because of the differential impacts resulting from different upstream watershed conditions that send FIO loads downstream to the coast.
In addition to the enteric pathogens of human health concern that can be indicated using FIO, microorganisms that cause other health problems, such as respiratory and skin infections, are also often present at beaches, and may be tracked to upstream sources [131]. This is especially true in the context of extreme storm events when beaches are heavily impacted by upstream river flows and plumes. Methods such as microbial source tracking (MST) and the monitoring of eDNA have shown value in their ability to improve predictive modeling of extreme storm events by offering insights into sources and transport pathways for FIO [7,132]. Differences in MST and eDNA monitoring results between "normal" and heavy storm conditions can be helpful in determining the types of microorganisms that become active within the aquatic environment in response to storm conditions [133]. Similarly, they can be informative in characterizing upstream impacts on coastal areas, by revealing potential catchment sources of microorganisms. Integrating well-calibrated watershed FIO models with nearshore water quality models (e.g., Bedri et al. [53]) or statistical and data-based approaches that describe FIO loading to coastal areas [20] may further improve the performance of nearshore FIO models during extreme events.
Conditions surrounding FIO sedimentation, attachment to suspended solids, and resuspension in riverbeds and coastal areas can vary greatly between storm events. However, there is a notable lack of observational data on water quality during and immediately following different storm conditions. High-resolution FIO data both within and between storm events will be critical to effective simulation of FIO loading, attachment dynamics, sedimentation and resuspension kinetics, and overall water quality in river plumes associated with heavy rain events. It can be difficult to collect these data, due to safety considerations, but the use of sensor networks and small unmanned aerial vehicles has emerged as a potential alternative to field data collection. Morgan et al. [134] demonstrated the use of unmanned aerial vehicles to photograph and document inland irrigation ponds and used image analysis to characterize water quality from the images collected. Many existing sensors on water bodies (e.g., select USGS gauging stations) also collect ancillary water quality data such as turbidity and electrical conductivity. These easy-to-collect data have the potential to help further constrain and evaluate FIO models because of their correlations to microbial water quality [10,11,16,135]. High-resolution water quality data collection is ideal for effective fate and transport modeling, but this data collection can take many forms, including remote sensing and proxy data collection.
In the coastal environment, there are multiple ways to simulate river plumes in numerical water quality models. Along with the river flow inputs, plumes may be characterized by tracking specific FIO within a water quality or FIO-specific model. In these models, FIO concentrations associated with the river inputs and decay function parameters can be specified to reflect local conditions. To assess the relative contributions of FIO from riverine sources to a beach site, constant FIO concentrations or arbitrary FIO masses can be input into the model over a release period [136,137] and breakthrough curves can be generated over time for specific locations. In contrast, FIO concentrations that are associated with riverine flows may be calculated using empirical relationships between river flowrate and FIO concentration [16,17,20]. For beaches impacted by multiple river plumes (e.g., Liu et al. [13]; Kim et al. [138]), the plume dynamics and hence nearshore water quality can be significantly more complex (Figure 1). Using realistic boundary conditions/forcing, models can track the FIO within the plumes spatiotemporally. Another attractive option for simulating FIO plumes involves the use of particle tracking [96,[139][140][141][142], especially "reactive" particle tracking models that can account for FIO losses [143]. In this case, FIO are released to the model domain (i.e., river outlets) as discrete particles. Upon their release, the particles' movements are tracked over time based on the simulated velocity field in three dimensions. This approach, using a Lagrangian formulation for dispersion, has the advantage that it does not suffer from excessive numerical dispersion inherent to Eulerian approaches. All of these approaches have their merits and drawbacks, so it is likely that selection of an optimal framework for plume modeling will require evaluation of the approaches within the context of the research questions and local conditions. Emerging issues such as water quality degradation associated with FIO exchange between water and sand, river plumes, upstream watershed impacts, and heavy storm runoff are key to effective modeling of microbial fate and transport in coastal areas. As such, future research and modeling in these areas will be beneficial to the water quality modeling community and knowledge base into the future.

Conclusions
Numerical models of water quality in nearshore regions can be useful tools for management of recreational water resources. Water quality and FIO fate and transport models within larger coastal ocean modeling frameworks have the potential to predict the fate and transport of FIO and pathogens of human health concern over time and across space. However, these models are only useful if they are refined and validated against observations.
In recent years, development of reliable FIO fate and transport models for aquatic and coastal systems has been an active area of research. As a result, modeling approaches frequently include decay terms associated with base mortality, solar inactivation, and sedimentation, though the specific parameterization of those terms can vary between models [41]. Highly generalized fate and transport models expand upon those terms to account for the effects of salinity, water temperature, pH, dissolved organic carbon, and differential settling rates due to varying particle sizes [21]. Model optimization in terms of FIO tracking has led to frameworks with RMSE values as low as 0.41 log 10 FIO CFU 100 mL −1 of water [15]. Some models have also been shown to predict up to 87% of variation in FIO concentrations from observed data [18]. While these validation statistics indicate that model frameworks are improving in their prediction of water quality, there is still room to optimize further. It is also important to note that many of these model parameterizations are specific to their local model domains. Generic models of FIO fate and transport can be developed, but without extensive datasets to test and constrain processes generic model formulations may not offer superior performance compared to simpler models [21]. Therefore, it will likely continue to be imperative that models be developed for their specific contexts, in order to maximize their predictive capacity.
FIO fate and transport modeling frameworks linked to watershed models in the contexts of water-sand exchange at the beach and release of FIO during storms can help us prepare for the potential impacts of extreme events on coastal areas. High quality intra-and inter-event data as well as modeling studies are needed to push the predictive capability of the current generation of FIO models. By refining established FIO decay functions to maximize predictive ability of models and combining those with the diffuse point and non-point FIO sources like plumes and sand-water exchange, prediction and tracking of pollutants in nearshore water and sand can be optimized. Confidence in modeling results can be maximized, allowing for more effective management for public health at nearshore and recreational beach areas in the face of climate and land use change.