Progress and Challenges in Coupled Hydrodynamic-Ecological Estuarine Modeling

Numerical modeling has emerged over the last several decades as a widely accepted tool for investigations in environmental sciences. In estuarine research, hydrodynamic and ecological models have moved along parallel tracks with regard to complexity, refinement, computational power, and incorporation of uncertainty. Coupled hydrodynamic-ecological models have been used to assess ecosystem processes and interactions, simulate future scenarios, and evaluate remedial actions in response to eutrophication, habitat loss, and freshwater diversion. The need to couple hydrodynamic and ecological models to address research and management questions is clear because dynamic feedbacks between biotic and physical processes are critical interactions within ecosystems. In this review, we present historical and modern perspectives on estuarine hydrodynamic and ecological modeling, consider model limitations, and address aspects of model linkage, skill assessment, and complexity. We discuss the balance between spatial and temporal resolution and present examples using different spatiotemporal scales. Finally, we recommend future lines of inquiry, approaches to balance complexity and uncertainty, and model transparency and utility. It is idealistic to think we can pursue a “theory of everything” for estuarine models, but recent advances suggest that models for both scientific investigations and management applications will continue to improve in terms of realism, precision, and accuracy.


Introduction to Models
Models are defined in numerous ways: Ba usually small copy of something,^Ba particular type or version of a product,^and most relevant to environmental models, Ba set of ideas and numbers that describe the past, present, or future state of something^ (Merriam-Webster 2014). While the latter most embodies simulation models, the other definitions allude to the role of models as versions of real (or hypothetical) systems. Models can involve ideas or numbers that describe a system, thus encompassing both conceptual and mathematical constructs. Reckhow and Chapra (1983b) define a model as Ba simplified representation of a real object, process, concept, or system^while Haefner (2005) defines them generally as Ba description of a system.^Many authors have coined their own definitions, but all highlight models as idealized representations of natural systems. Brush and Harris (in press) define models as Bconceptual or mathematical simplifications of real systems.^This includes conceptual models, in which relationships between system components are diagrammed to synthesize understanding of system function and enable qualitative predictions of future Communicated by Iris C. Anderson states (e.g., Valiela et al. 1997a;Jay et al. 2000;Thomas et al. 2006). Mathematical models may explore qualitative system behaviors and have been critical in synthesizing data and decision making in ecologically sensitive regions (Di Gennaro et al. 2012). Peters (1991) emphasized development of falsifiable, quantitative predictions in ecology, where models are used as tools to test logic based upon a defined set of assumptions. Models have moved beyond the conceptual phase and involve predictive equations to describe real phenomena; these models can take the form of empirical or statistical relationships driven by observed data (e.g., Dillon and Rigler 1974;Vollenweider 1976;Nixon et al. 1996Nixon et al. , 2001 or mechanistic, process-based models that piece together equations representing underlying physical, chemical, and biological processes (e.g., Kremer and Nixon 1978;Cerco and Cole 1994;Cerco and Noel 2004;Luff and Moll 2004;Soetaert and Middelburg 2009). Some mechanistic models have exact, analytical solutions (e.g., Streeter and Phelps 1925;Hansen and Rattray 1965), but most must be discretized in space and time and solved by numerical simulation. It is the latter type of model that is the focus of this review.
In both estuarine hydrodynamics and ecology, the use of simulation models has increased during the past several decades Harris 2010, 2015). Numerical models of estuarine hydrodynamics have expanded from tidally averaged onedimensional models (Hansen and Rattray 1965) to fully threedimensional (3D) models that can resolve wetting and drying, wave-current interaction, and sediment transport (e.g., Warner et al. 2005;Zhao et al. 2010;Olabarrieta et al. 2011). Numerical models emerged as heuristic tools for understanding ecosystem structure and function in the 1970s (e.g., Patten 1971;Di Toro et al. 1971;Odum 1971;Kremer and Nixon 1978) and their utility as management tools grew in the 1980s (e.g., HydroQual 1987). The development of increasingly powerful computers, user-friendly software, and the acceptance of models as research tools led to expanded use in the 1990s and the current century. With widespread use in both the research and management arenas, a review of the history, applications, limitations, and future directions of estuarine models is timely.

Historical Trajectory of Estuarine Hydrodynamic Models
Numerical solutions to estuarine hydrodynamics have existed since the pioneering work of Hansen and Rattray (1965). At the simplest level, one-dimensional vertical models solve the Reynolds-averaged Navier-Stokes equations for fluid flow, with numerous idealizations, and yield a vertical profile of velocity and salinity. These 1D modeling exercises recreated the net landward transport of water near the bed and net seaward transport at the surface, i.e., estuarine circulation (see Geyer and MacCready 2014 for a comprehensive discussion). Practically, these solutions had important ramifications for understanding the transport of salt, sediment, and contaminants.
Inclusion of the longitudinal (along-estuary) coordinate extended the models into two dimensions and allowed for the identification of zones of estuarine convergence and stagnation. Festa and Hansen (1976) used such a model to identify how salinity intrusions respond to changes in depth and river flow, demonstrating a landward movement of an intrusion with channel deepening. While this treatment was sufficient for narrow and uniform estuaries, systems with an irregular geomorphic layout required inclusion of the lateral dimension, leading to a class of two-dimensional, depthaveraged models. Cheng et al. (1993) extended the depthaveraged approach in a model of San Francisco Bay to include wetting and drying of intertidal areas. This model was instrumental in understanding spatial patterns in tidal velocities and residual circulation that previously could not be resolved by observations alone. Officer (1980) developed an alternative approach for simulating water, salt, and constituent transport through the use of box models. These models use forced freshwater inflows and salinity distributions to solve for transport in a series of homogenous boxes along the estuary axis through a salt balance; boxes can be depth-averaged or stratified (two layers). In onelayer systems, Officer box models simulate the seaward flow of water; in two-layer systems, they represent the gravitational circulation, predicting seaward flow in the surface layer, landward flow in the bottom layer, and vertical entrainment of salt between the layers. Box models can also include dispersive tidal transports using the Hansen-Rattray parameter to specify tidal dispersion (Officer and Kester 1991).
Ultimately, many estuaries are not conducive to lateral-or depth-averaging, and 3D solutions may be required for accurate representation of hydrodynamics due to complex bathymetry, vertical gradients, or lateral circulation. For example, Warner et al. (2005) developed a 3D model of the Hudson River estuary to investigate variations in stratification during spring and neap tides, demonstrating the skill of the model in capturing tidal and subtidal timescale processes. Ralston et al. (2010), simulating the Merrimack River estuary using an unstructured grid model, showed that resolution of bathymetric variability has a large effect on proper simulation of salt fluxes and stratification.
Increases in computational resources have allowed for development of more robust mixing routines, higher-resolution domains, and inclusion of detailed circulation processes. Early models relied on specified vertical mixing, but advanced computing power has allowed for two-equation turbulence models that solve for the transport of turbulent kinetic energy and turbulent length scales (Umlauf and Burchard 2003). Advanced horizontal transport routines such as MPDATA (Smolarkiewicz and Margolin 1998), which are critical for robust simulation of scalars and constituents, are sizeable improvements over earlier advection routines. In terms of resolution, Zhao et al. (2010) implemented an unstructured grid model of the Plum Island estuary with a maximum horizontal resolution of 7 m in tidal wetland creeks to investigate inundation processes in marshes while also simulating salt transport on the inner shelf. Increased computational power can also uncover novel hydrodynamic mechanisms: Olabarrieta et al. (2011) modeled wave-current interaction with the vortex force method in Willapa Bay, Washington and found that wave-induced setup during storms affected exchange between the estuary and coastal ocean.
The increase in computational capacity is clear in the exponential increase in computational elements and grid resolution across numerous modeling systems (e.g., ROMS, Delft3D, SELFE, FVCOM, UnTRIM, TRIM3D, EFDC), computational schemes (finite-difference, finite-element, finite-volume), and domain discretizations (rectilinear, curvilinear, unstructured, nested) (Fig. 1). This trend of increasing resolution reflects the realization that complex bathymetry and geomorphic layout have a pronounced effect on circulation (and thus biology) and should be represented. Given the availability of computational resources, most researchers therefore opt for increased resolution wherever possible. As the following section demonstrates, estuarine ecological modeling has followed similar trends, though there is also a realization that reduced complexity may result in a more robust simulation.

Historical Trajectory of Ecological Modeling
Ecological models of aquatic systems can be traced back to the pioneering work of Gordon A. Riley (1946Riley ( , 1947. Riley adapted the Lotka-Volterra predator-prey equations to develop the first mechanistic, process-based models of phytoplankton (Riley 1946) and zooplankton (Riley 1947) in a marine system, using Georges Bank in the North Atlantic as a test site. Riley's major advancement was to formulate models in which changes in daily growth and loss processes were functions of biomass, abundance of prey and predators, and environmental variables such as temperature, light, and nutrients. These independent models of phytoplankton (P) and zooplankton (Z) were then combined along with a state variable for nutrients (N) to produce the first NPZ model of an aquatic system (Riley et al. 1949).
The next major phase in the application of ecological models to aquatic systems came during the 1960s and 1970s with the development of models primarily for estuaries, shelf ecosystems, and lakes (summarized in Brush and Harris 2010;Fig. 2). Early examples included Di Toro et al.'s (1971) model of the Sacramento-San Joaquin Delta, Steele's (1974) model of the North Sea, and Kremer and Nixon's (1978) model of Narragansett Bay. These early models had limited biological and spatial resolution, i.e., one to two phytoplankton groups, one zooplankton group, cycling of a few inorganic nutrients, and were either zero-or one-dimensional with relatively coarse spatial elements. The use and complexity of models continued to expand in the 1980s to systems such as Chesapeake Bay (HydroQual 1987), the Baltic Sea (Stigebrandt and Wulff 1987), and the Ems-Dollard estuary (Baretta and Ruardij 1988). These models added complexity, with dissolved oxygen as a state variable and addition of dissolved organic and particulate forms of nutrients and implementation in higher-resolution, 3D domains. The Ems-Dollard model included pelagic, benthic, and epibenthic food webs, including microphytobenthos, meiobenthos, deposit and suspension feeders, and two groups of epibenthos. The application of models to management began in this period, particularly related to anthropogenic nutrient enrichment and the cultural eutrophication of estuaries (HydroQual 1987).
Models have evolved with an increasing focus on management (HydroQual 1991;HydroQual and Normandeau Associates 1995), increasing hydrodynamic resolution (Fig. 1), and addition of biological and biogeochemical detail, including watershed-estuary interactions and numerous components of coastal ecosystems and their food webs (e.g., Cerco and Cole 1994;Baretta-Bekker 1995;Baretta-Bekker and Baretta 1997;Billen and Garnier 1997). For Chesapeake Bay, Cerco and Cole's (1994) model included 22 water column state variables including partitioning of organic material into labile and refractory pools and was coupled to a detailed sediment flux model with ten state variables (Di Toro 2001). The 2002 version of the Chesapeake Bay model operated in nearly 13,000 cells (∼ 2-km horizontal resolution; Cerco and Noel 2004), which has since been increased to 57,000 (US EPA 2010); it also includes three species of seagrass, micro-and mesozooplankton, deposit and suspension feeders, Eastern  Fig. 1 Trends in estuarine hydrodynamic modeling complexity, from a computational cell and grid resolution perspective. Models were selected from peer-reviewed publications that contained full reporting of grid domain characteristics. A variety of models are represented spanning several continents, modeling systems, and estuarine types oysters, and Atlantic menhaden (Cerco and Moore 2001;Noel 2004, 2007;Dalyander and Cerco 2010). Models are now essential management tools (e.g., US EPA 1999;NRC 2000;Giblin and Vallino 2003;Harris et al. 2003;Wulff et al. 2007) and are used to develop total maximum daily loads (TMDLs) for numerous estuarine pollutants (e.g., US EPA 2010).
Despite the trend towards increasing model complexity and resolution, there are trade-offs among complexity, resolution, over-parameterization, and model uncertainty. Numerous investigators have addressed the role of model complexity (Fulton et al. 2003(Fulton et al. , 2004Ménesguen et al. 2007) and others have embraced the use of reduced complexity and alternative modeling approaches (Rigler and Peters 1995;NRC 2000;Pace 2001;Duarte et al. 2003;Scavia et al. 2006;Swaney et al. 2008). Multiple modeling approaches have also been combined to inform coastal management (e.g., Stow et al. 2003;Scavia et al. 2004). Clark et al. (2001) promoted abandoning prediction as a goal for ecological models, instead focusing on forecasting as a more pragmatic modeling objective. However, this viewpoint reverts to a pessimism that was warned against by Peters (1991), who advocated for quantitative prediction as the hallmark of a robust ecological model. Harris et al. (2006) proposed leveraging the metabolic theory of ecology (MTE) to understand and predict estuarine ecosystem function. The MTE is notable in ecology because predictions and derivations may arise from first principles alone providing a parallel to the Newtonian physics used to describe estuarine hydrodynamics. Any efforts in ecology that improve the theoretical basis of our models and predictions will be critical to forging models that combine these two fields. Harte (2002) acknowledges the challenges of combining these scientific perspectives, calling on a synthesis of Newtonian and Darwinian traditions. Recent examples to leverage the MTE in ecological models (Cheung et al. 2008;López-Urrutia 2008;Harris and Brush 2012;Sinsabaugh et al. 2013;Ehrich and Harris 2015) hint at the potential for accomplishing this goal in estuarine models.

Philosophical Considerations in Numerical Modeling
Modelers and users of model-derived knowledge are burdened by the expectation that model results should mimic the world, while also advancing and aligning with theory. Indeed, theoretical physicists pursuing the standard model of elementary particle behavior allude to defining a Btheory of everything.^In practice, those working at the ecosystem scale have described the constraints of numerical models designed to explore ecology (e.g., Levins 1966), hydrodynamics (e.g., Blocken and Gualtieri 2012;Roache 1997), and biogeochemistry (e.g., Anderson 2005). Morrison and Morgan (1999) identify three features of models that clarify their role: (1) they are constructed in partial independence and dependence on both the world and theory, (2) they function autonomously from the world and theory, and (3) they are capable of representing aspects of the world and theory at the same time. Similar to the idea of conformation, where a representation of a system conforms to some aspect of the system (i.e., a road map of Albuquerque can reliably represent the roads of Albuquerque, but not the distribution of lizard habitat in Albuquerque), reliable models can only be expected to conform to the discretized equations they are solving (theory) or type of observational data to which they are calibrated (reality + measurement noise). A reliable model cannot be expected to represent every aspect of the observational data (e.g., temporal resolution of the model is coarser than real time) nor is the representation exactly accurate as there is always some quantifiable error. Longino (2002) writes, BConformation is more suitable than true or false for expressing the ways in which Year Fig. 2 Trajectory of aquatic ecological modeling as illustrated by the number of publications using the term Becosystem model^in the Aquatic Sciences and Fisheries Abstracts (ASFA) database. Modified from Harris (2010, 2015) complex content, such as a theory or model, is successful representation.^Therefore, a model that reliably conforms to some aspect of the system of interest can be considered a successful and reliable representation of that aspect. Again, reliability indicates that the model consistently and accurately (within some error tolerance) represents the equations being solved and the observations of interest. Levins (1966) outlined the competing goals of generality, precision, and realism. Modelers may strive to make their models general, in which they are applicable to a variety of systems with little adjustment; precise, in which they reproduce the observed data with minimal error; and realistic, in which they attempt to explain the dynamics of the real system in detail. Unfortunately, these three goals cannot be simultaneously achieved in practice. Increasingly, realistic models would include as many processes and components as possible; however, the inclusion of poorly constrained parameters likely erodes model precision. Simple empirical/statistical models may maximize precision (Rigler and Peters 1995); however, they lack the processes providing explanatory power and cannot be used to extrapolate beyond the known range of values. General models may need to be optimized in a way that sacrifices some level of both precision and realism to be broadly applicable. Such decisions are often driven by management context or modeling goals. Formal methods exist for identifying optimal levels of model complexity for single formulations or statistical models based explicitly on the information in the data (Burnham and Anderson 2002), Similar approaches for multi-parameter, multi-state variable models are beginning to emerge in the Bayesian arena (e.g., Obenour et al. 2014).

Parsimony and Over-parameterization
Efforts to maximize detail and utilize the highest degree of computational power often lead to excessive parameterization, greater model uncertainty, and error propagation (Denman 2003;Anderson 2005;). On the other hand, simple models with few processes can lack important feedbacks and interactions that are necessary for management applications (Reckhow 1999;Flynn 2005;Evans et al. 2013). To overcome these challenges, it is important to first delineate the conceptual extent of the model by focusing on (1) the research question one intends to address, (2) the current understanding of critical processes and the study system, (3) the availability of data to support parameterization, and (4) how the model will be used in the future.
Ideally, the research question should define the level of realism required of a model, which may range from simple empirical formulations to coupled models based on differential equations. Often, however, processes and state variables beyond the scope of the research question are included to explore other ecological processes. Adding processes and variables increases the number of parameters and can quickly surpass our ability to constrain the formulations (Denman 2003). If, for instance, our research question is focused on predicting primary production, this might be achievable using only five variables in an empirical model, as opposed to a fully coupled atmospheric-watershed-estuarine model (see Fig. 3). It is also important-though sometimes overlooked-to ensure that there is an adequate amount of data to calibrate and validate the processes represented in the model and any related formulations (Flynn 2005).
Several inter-comparative modeling studies and reviews have considered the degree of complexity necessary to best represent fundamental biogeochemical processes Raick et al. 2006;Friedrichs et al. 2007) and multi-species trophic interactions (Fulton et al. 2003). These comparisons frequently focus on multiple criteria including performance-and skill-based metrics (e.g., correlation, RMS difference, standard deviation, and cost functions). While the degree of model complexity in these studies ranged from relatively simple NPZ models to complex ecosystem models (containing multiple phytoplankton, zooplankton, nutrient, and detrital pools and associated  Fig. 3 The development of ecosystem models requires a balance between inclusion of critical processes while limiting the addition of unconstrained formulations. For example, a model focused on predicting primary production may only require the five key variables located closest to the inner circle. Including additional complexity, like the processes located in the outer circles (e.g., dynamic nutrient cycles, multiple consumer groups) may account for other mechanisms we wish to capture. However, their inclusion will introduce new formulations into the model which may not be well constrained. CDOM chromophoric dissolved organic matter, DIC dissolved inorganic carbon, DIN dissolved inorganic nitrogen; DIP dissolved inorganic phosphorus, DOC dissolved organic carbon, POC particulate organic carbon; K D : light attenuation, TSS total suspended solids biogeochemical cycling), the results highlighted that simple and reduced complexity models were able to reproduce data for specific regions for which they were calibrated with the same skill as more complex models (Raick et al. 2006;Friedrichs et al. 2007). The broad use of models for management and policy applications and the limitations outlined above highlight the value of reduced complexity models (Pace 2001;Duarte et al. 2003;Ménesguen et al. 2007, Lucas andThompson 2012). The need for balancing parsimony with over-parameterization leads to three commonly observed approaches: (1) begin with a biologically and physically reasonable simple model, and add complexity as needed; (2) begin with a complex model and remove complexity; or (3) use multiple models with varying levels of complexity, which can be used as ensembles to predict a range of future states (e.g., Stow et al. 2003;Scavia et al. 2004). The latter approach recognizes that no single model can comprehensively represent the processes and goals for all potential applications.
This review, though focused on the future of estuarine hydrodynamic and ecological modeling, describes a diversity of models. In the case of numerical models developed to meet the needs of management for TDMLs of a pollutant to an estuary, many use systems of ordinary differential equations that describe the assumed mechanisms most pertinent to the desired water quality restoration goals. These models commonly use the pollutant (such as nitrogen or phosphorus) as a currency to drive these criteria so that a quantitative loading goal may be predicted. These management tools contrast with models developed to explore the system dynamics governing ecosystem feedbacks or regime shifts. Examples of the latter include models that successfully predict regime shifts of submerged aquatic vegetation (SAV; Carr et al. 2012;del Barrio et al. 2014) and shellfish larval recruitment success (Bidegain et al. 2013). A theory of everything for estuaries may be neither practical nor desirable, but recent modeling advances suggest that modeling tools for either science or management will continue to improve in realism, precision, and accuracy. In the meantime, we must be clear with managers and other end-users regarding the appropriate application and inherent limitations of estuarine hydrodynamic and ecological models.

Model Limitations and Their Influence on Simulating Ecological Processes Hydrodynamics
For numerical models, the equations governing estuarine hydrodynamics are well constrained and tractable when compared with ecological processes. Correspondingly, field measurements of physical quantities such as water level, water velocity, temperature, salinity, and PAR are more easily acquired with modern instrumentation than ecological variables such as nutrient or organism concentrations, feeding rates, or biogeochemical transformations. Nonetheless, data on physical quantities are not always available at the necessary spatial scales. Ecological model results are quite sensitive to the hydrodynamic framework Allen et al. 2007) through advection and dependence on properties such as temperature or irradiance, so robust evaluation of the hydrodynamic component is imperative for coupled models.

Spatial Resolution
Due to computational limitations, a trade-off exists between spatial resolution and simulation duration. Correspondingly, the ability of models to represent physical processes depends on the spatial resolution. A hydrodynamic model needs to be able to represent the scales of physical processes relevant to the ecological model, particularly for non-linear processes that are not well represented by spatial or temporal averaging. Strong gradients in hydrodynamic properties, such as salinity fronts or those associated with lateral circulation, can correspond with ecological patchiness that would not be represented in lower-resolution models. Simulation over timescales from years to centuries requires compromise in spatial resolution (Penduff et al. 2010) and parameterization of sub-gridscale hydrodynamic and ecological processes.

Turbulent Mixing
Few hydrodynamic and ecological models resolve the spatial and temporal scales of turbulence; thus, some form of parameterization is needed to represent mixing at scales smaller than the grid size. Turbulent mixing can be represented simply with constant eddy diffusivity or with a two-equation turbulence model (e.g., Mellor-Yamada, k-epsilon, or k-omega). The degree of complexity (constant, one-equation, two-equation) of the mixing scheme has greater consequences for the ecological model results than variations among the different twoequation turbulence closure schemes (Umlauf and Burchard 2003;Warner et al. 2005). In most estuaries, stratification due to salinity or temperature significantly modifies rates of vertical mixing for at least part of the tidal or seasonal cycle, and thus, it is necessary to incorporate the effects of stratification on transport.
Turbulence closures impose rates of vertical mixing on scalars, but additional mixing occurs from numerical diffusion due to advection of scalar gradients. Rates of numerical mixing depend on the advection scheme, on the grid resolution, and on the scalar gradients (Burchard and Rennau 2008). In models with lower order advection schemes, coarse discretization, or strong scalar gradients, the hidden numerical mixing can be similar to or greater than the physical mixing from the turbulence closure. While a minimum or background diffusivity is often imposed to reduce instabilities, numerical mixing can be the limiting factor in the ability of models to represent sharp gradients such as pycnoclines, lutoclines, or thin phytoplankton layers.
Many estuaries are relatively shallow, and turbulence is generated predominately by shear at the bed. The gravitational exchange flow due to the salinity gradient is also extremely sensitive to water depth. Consequently, accurate representation of the bathymetry is critical for hydrodynamic models (Ralston et al. 2010, Ganju et al. 2011. In coarser resolution models, the projection of bathymetry onto the grid can affect the hydrodynamic results. For example, in a cross-sectionally averaged model of circulation and salinity in the Hudson River estuary, the channel thalweg proved to be the relevant depth for predicting the salinity distribution, but modeling sediment transport required both the channel and shoal bathymetry (Ralston and Geyer 2009). At any model resolution, highquality bathymetric data is key to reducing uncertainty in the hydrodynamic model (Plant et al. 2009).

Bottom Shear Stress and Sediment Transport
The distribution of bottom shear stress depends on bathymetry, water velocities, and bottom roughness. Parameterization of the bottom roughness, whether with a roughness length (z 0 ) or drag coefficient (C d ), is a major source of uncertainty for the hydrodynamic model (Davies and Lawrence 1995;Ganju and Sherwood 2010). Bottom roughness typically depends on bed composition, but rarely are maps of grain size, bed forms, or vegetation available at the model grid resolution. Instead, bottom roughness is often made uniform, defined based on empirical relationships, or made spatially variable as a tuning parameter. None of these options is robust, and the associated uncertainty affects vertical mixing, sediment erosion and deposition, and benthic fluxes such as dissolved oxygen or larvae. Uncertainty in the bottom roughness can be reduced through bottom stress measurements, for example, with acoustic Doppler velocimeters, but collection and analysis of such measurements are difficult.
In addition to the dependence on bottom roughness, modeling estuarine sediment transport is highly dependent on the sediment properties (Amoudry and Souza 2011), which are typically underconstrained by observational data. Important sediment model inputs include the number of size classes, settling velocities, erosion rates, critical shear stresses, and initial bed conditions. Model setup and evaluation requires characterization of sediment properties both in the bed and in the water column, and sediment properties may vary temporally with storm events or biological activity. Fluxes of sediment between the bed and water column can directly affect ecological model results through light attenuation De Boer 2007) and other water-column processes.

Specification of Lateral Boundary Conditions
Defining hydrodynamic boundary conditions in most estuarine models is relatively straightforward compared with the open ocean and often includes tidal water level and river discharges. However, boundary conditions can be constrained by data availability. In estuaries where wind affects mean flows and mixing, mainland weather stations may inadequately represent the wind magnitude and distribution over water (Scully 2013), and an atmospheric model that resolves the wind distribution at the scale of the estuary may be necessary (Raubenheimer et al. 2013). In estuaries without a major river, groundwater may provide the primary source of freshwater and nutrients, but determining groundwater fluxes requires major observational effort or alternative modeling analysis (Ganju et al. 2012).

Ecology
Biological systems are inherently complex and several processes are still poorly understood. Our capability to model these processes is still notably limited by observational data availability as well as computational expense. These limitations illustrate the difficulty in balancing generality, precision, and realism; they also show how parsimony and overparameterization compete in ecological modeling applications.

Complexity of Biological Systems
Compared to hydrodynamic models, ecological models are limited by the lack of fundamental deterministic equations. As Harte (2002) stated in his comparison of the physical and ecological worldviews, BPhysicists seek simplicity in universal laws. Ecologists revel in complex interdependencies. A sustainable future for our planet will probably require a look at life from both sides.^Our conceptual understanding of marine food webs has evolved from simple food chains to complex, interconnected webs with numerous species, life histories, and generalists that feed at multiple trophic levels (Landry 1977;Baird and Ulanowicz 1989;Nixon 1992). A model cannot include all species or functional groups, which therefore necessitates aggregation into larger groups and the inclusion of closure terms (i.e., specification of grazing by one trophic level higher than that which is modeled) to ensure model stability . While the tendency has been towards increasing complexity through inclusion of more state variables and processes, our understanding of many of these variables is still incomplete, as are the data to parameterize and verify the models. Ecological systems are characterized by high levels of individual variation, weak explanatory relationships, strong dependence on previous events, and numerous potential future states that are difficult to predict (Low-Décarie et al. 2014).

Biological/Biogeochemical Process Understanding
Several factors are recognized as gaps in our understanding of water quality modeling in estuaries. In a recent review, Statham (2012) noted that improved understanding of cycling of organic N and P, resuspension, benthic exchanges, submarine groundwater discharges, and rates and magnitudes of bacterially driven processes is needed. Challenges for endto-end marine modeling, which simulates systems through higher trophic levels, include representing zooplankton as a link between lower and higher trophic levels; simulating macroinvertebrates, demersal fishes, and behavioral movement of fishes; and including feedbacks from biota to physicalchemical processes (Rose et al. 2010). To address these gaps, it is important to develop methods and studies that can advance our understanding of poorly constrained model formulations. Uncertainty is associated with model structure and complexity; Link et al. (2012) provided multiple methods to address this. In addition to these recognized aspects of uncertainty, the model may contain hidden uncertainties which may not be fully characterized. The iterative nature of adaptive management allows for system models to incorporate additional information and understanding as it becomes available (Stow et al. 2003), providing a means of addressing uncertainties as they are recognized. For many processes, understanding of basic ecological dynamics has not advanced sufficiently to be accurately incorporated into models (e.g., Anammox; van Niftrik and Jetten 2012). Given the difficulty of sampling processes beyond limited spatiotemporal scales, a relevant question arises: should we ignore poorly quantified processes or use sparse data to parameterize the best we can?

Spatial and Temporal Resolution
The spatial scales of processes often drive the spatial resolution in a model. In some estuaries, regional patterns tend to dominate seasonal biogeochemical cycles. Much of the data available to validate models are collected at regionally based monitoring stations or measured intensively at a few strategic locations during short-term research studies. Thus, there is a utility in using regional scale model computations (i.e., box and tidal prism models) in many ecosystems (Lake and Brush 2015) and when developing and calibrating ecological formulations (Ménesguen et al. 2007). However, continuous monitoring by buoy arrays and other observation networks has demonstrated patchiness in phytoplankton production and biomass (Harding et al. 2005;Martin 2003). Patchiness can exist in perfectly uniform physical environments through mechanisms theoretically attributable to diffusive instabilities (e.g., Mimura and Murray 1978) or social interactions (e.g., swarming, Okubo 1980). Although this suggests that physical structure is not necessarily predictive of biological structure at all scales of interest, modelers have sought to increase spatial resolution to capture dynamics at these finer scales. It has been argued, however, that Bsmall-scale temporal and spatial resolution gives the illusion of substantial knowledge^ (Reckhow 1999), when we simulate spatial patterns that we cannot confirm, cannot validate, and may not fully understand.
Temporal resolution is typically less limiting for biogeochemical processes than spatial resolution. The timescale of most highly resolved models (e.g., minutes) are sufficient to capture the dynamics of the key biogeochemical processes. In fact, even the earliest marine ecosystem modeling studies were able to capture seasonal cycles of phytoplankton and zooplankton dynamics (Riley 1946). In many instances, the underlying formulations represent a daily approximation (Kremer et al. 2010) and are based on experiments measured on timescales of a few hours. This raises several important considerations: (1) is it realistic to scale measured rates to such short model time steps (minutes), (2) are model time steps limited by past measurements and sampling techniques, and (3) what previous model formulations (e.g., phytoplankton growth rates) may incorrectly predict biomass when applied at short time steps (Ménesguen et al. 2007)? In many instances, high temporal resolution is required to capture short-term processes including dynamic chlorophyll and oxygen cycles in shallow, productive systems where concentrations can vary by orders of magnitude in a given day (e.g., Tyler et al. 2009) or in systems where aggregation of hydrodynamic output (to a 1-h time step) may fail to resolve tidal dynamics.

Linking Hydrodynamic and Ecological Models
Estuarine hydrodynamics and ecology can be simulated using combinations of disassociated or integrated models and may be run simultaneously or sequentially. We must distinguish first between these modeling configurations and their benefits and drawbacks. The first type of configuration is the offline, no-feedback modeling system (Fig. 4). In this case, hydrodynamic simulations are run for a prescribed period, the output is stored, and then it is used as input to a separate ecology model with or without spatial or temporal modification (e.g., Wool et al. 2003). One benefit is that the ecological model can be run multiple times with the same hydrodynamic input, eliminating the need for repeating costly hydrodynamic calculations. One drawback is that some post-processing of hydrodynamic model output may be necessary, and this may change the fidelity of the hydrodynamic results (e.g., vertical mixing may not be congruent with the original simulation). In the second configuration, online with no feedback, the models are integrated into the same compiled executable, such that hydrodynamic model output is passed internally at runtime to the ecological model with no need for modification. There is no loss of fidelity in this case, because the models are ostensibly running on the same grid and time step; however, one must re-run hydrodynamics each time an ecological parameter is changed, even though the ecological parameter change has no effect on hydrodynamics. The last configuration is the online system with feedback between hydrodynamics and ecology: all routines are integrated into the same compiled executable, with feedback from the ecological model to the hydrodynamic model. The models may also be developed as modules that are coupled in a framework at runtime (Bruggeman and Bolding 2014).

Offline, No-Feedback Models
Offline, no-feedback models consist of two independent components that are integrated such that the output of one model is used as input for the next model. The one-directional passing of information may omit feedbacks between models, but it allows users to more quickly conduct sensitivity analysis, and run many scenarios.
Early examples of offline modeling are found in the work of Cloern (1991) and Koseff et al. (1993). In both of those studies, a one-dimensional (vertical) model for phytoplankton growth was implemented, with hydrodynamics incorporated only through the eddy diffusivity (turbulent mixing) term. Cloern (1991) used a vertically uniform value, while Koseff et al. (1993) prescribed realistic vertical profiles of eddy diffusivity. These studies showed how phytoplankton blooms in estuaries can be modulated by turbulent mixing and stratification, and explained the tendency of blooms to occur during neap tides in San Francisco Bay when mixing was minimized. In contrast to deep, partially mixed estuaries, shallow eutrophic estuaries often harbor extensive benthic primary producers. Trancoso et al. (2005) spatially integrated hydrodynamic output and repeated tidal forcing to drive a water quality module that added macroalgae as a primary producer in a eutrophic lagoon. Simulating hydrodynamics and ecology simultaneously would have prohibited annualscale simulations, and the model showed that inclusion of macroalgae produced water quality results more consistent with field data than with phytoplankton alone. Kremer et al. (2010) similarly analyzed hydrodynamic model output to develop a daily exchange mixing matrix for simulating ecology in Narragansett Bay. The approach provided a cost-and timeefficient solution that was appropriate when the spatial and temporal scales needed for the ecology simulation was coarser than the scale used for hydrodynamic simulation. Cerco and Noel (2013) modeled decadal-scale eutrophication processes in Chesapeake Bay with an ecological module that received offline input from atmospheric deposition, watershed, and hydrodynamic models. Decoupling the ecological component from the physical models allowed for a 21-year simulation that would likely be onerous in a fully online simulation. With this approach, the authors showed that annual variability in hypoxic volume and water quality was strongly controlled by physical processes such as stratification.
A different approach to offline, one-way coupling is the use of particle tracking in which numerical drifters follow the flow from the hydrodynamic model but are subject to biological behavior (vertical migration, growth, reproduction) or geological processes. For instance, Culberson et al. (2004) explored the effects of vertical migration on the dispersal of larval fish in the San Francisco estuary. North et al. (2005North et al. ( , 2008 conducted similar studies for striped bass eggs and oyster larvae in the Chesapeake Bay. More complex biological behavior including reproduction, growth, and mortality has been  Fig. 4 Examples of three model coupling configurations for simulating hydrodynamics and ecology; variables h, T, and S represent water depth, temperature, and salinity, respectively coupled to a hydrodynamic model of the San Francisco estuary to characterize the population dynamics of delta smelt (Rose et al. 2013) following an individual-based modeling (IBM) approach. Similarly, Canal-Vergés et al. (2014) simulated drifting macroalgae as evolving particles and studied their effect on eelgrass recovery in Odense Fjord (Denmark).

Online, No-Feedback Models
Online, no-feedback models implement the hydrodynamic and ecological algorithms simultaneously, but information only passes from one model to the other (typically from hydrodynamic to ecological). In this approach, hydrodynamic variables modify ecological processes, but there is no feedback from ecological variables to the hydrodynamics. This approach is computationally more expensive than the offline method, but ensures seamless passing of quantities with no spatiotemporal resampling or interpolation. For example, Lucas et al. (1999a, b) integrated phytoplankton as a scalar quantity into the advection-diffusion routine of a depth-averaged hydrodynamic model to investigate the effect of tidal processes on primary productivity. With their online model, they showed the importance of horizontal transport on bloom formation, extending the offline work of Cloern (1991) and Koseff et al. (1993). Online coupling can be optimized by updating the biogeochemical fields less often than the hydrodynamics: Xu and Hood (2006) modeled the phytoplankton and light attenuation changes in Chesapeake Bay during dry and wet conditions with a biogeochemical time step that was 20 times larger than the hydrodynamic component. This computational savings allowed for annual simulations and comparison of the relative influence of freshwater flow on ecological dynamics. Banas et al. (2007) modified the tracer transport scheme in a 3D hydrodynamic model to add chlorophyll as a scalar in simulations of Willapa Bay. They showed that most of the oceanic chlorophyll advected into the estuary leaves without being consumed by local grazing and tidal current structure.

Online, Feedback Models
Online modeling with feedback is necessary when ecological state variables affect hydrodynamics. One example is the influence of submerged and emergent vegetation on drag and circulation. In the case of emergent vegetation, increased drag traps sediment, increases accretion, and helps maintain elevation of the marsh plain relative to the water surface. Without this feedback in a model, deposition would likely be underestimated. SAV also tends to trap sediment, reduce turbidity, and increase light availability to seagrass beds. If this positive feedback is ignored in a model, light attenuation would be overestimated and productivity underestimated. Temmerman et al. (2005) implemented feedback between marsh vegetation and hydrodynamics in a 3D model of the Scheldt estuary. Drag and turbulence were modified as a function of stem diameter and number of stems per unit area. They found that this feedback led to the formation of tidal channels and the familiar levee-basin topography found in wetland complexes. Similarly, Chen et al. (2007) modified a 3D model to account for drag and wave attenuation due to seagrass. Reductions in shear stress and wave energy decreased sediment concentrations and transport, representing a positive feedback for maintenance of seagrass meadows. Future studies of emergent and submerged vegetation should account for feedbacks between biomass production and physical drivers. While sedimentary and biological light attenuation and absorption have been considered in several models (e.g., del Barrio et al. 2014), the effect of absorption as a reduction of available heat to alter stratification in the water column remains a research topic that could be addressed with an online coupled system (Wetzel et al. 2006;Lengaigne et al. 2007).

Skill Assessment Skill Metrics
Model assessment requires both qualitative and quantitative evaluation (Fitzpatrick 2009). Although qualitative assessments help answer heuristic questions about how Bgood^a model is, quantitative skill metrics are necessary for comparing, developing, tuning, and improving predictive models. Modeling is an iterative process and quantitative skill metrics provide objective means of evaluating the iterations. Stow et al. (2009) andFitzpatrick (2009) summarized many of the main approaches for skill assessment. Numerous statistical approaches have been used to quantify model goodness of fit on a univariate basis, including (1) basic error metrics (e.g., absolute, relative, and root mean square error), (2) correlation and regression of observations and predictions, (3) parametric (e.g., t tests) and nonparametric tests (e.g., Wilcoxin test), (4) plots of model bias, cumulative error, and cumulative density functions, and (5) receiver operator characteristic (ROC) curves (Thomann 1982;Reckhow and Chapra 1983a;Parrish and Smith 1990;Reckhow et al. 1990;Mason and Graham 1999;Haefner 2005;Brown and Davis 2006;Fitzpatrick 2009;Sheng and Kim 2009;Stow et al. 2009). Taylor and Target diagrams provide additional visual methods for displaying multiple error metrics on a single plot ).
Model skill scores provide a means of quantifying model goodness of fit; these are a measure of the normalized error of a model's predictive capability. Model skills are usually assessed by using the misfit between model predictions and data, but it is important to note that this misfit is not the same as the model error (Stow et al. 2009). In addition to the misfit caused by model error, uncertainty in how well the data reflects the true state of a system also contributes to the discrepancy between model predictions and observations. The model-data misfit is usually normalized, and how this is done greatly affects the skill score and what it represents.
An instructive approach to skill assessment is the Brier skill score (Brier 1950), which normalizes the model misfit by the misfit relative to a Breference^model: where the vector of observations x o,i contains N measurements, x p,i are the model predictions at the same time and location of the data, and the vector x r,i are the predictions of the Breference model.^The reference model may be climatological data, or if the reference model is taken to be the mean of the data vector, x o , then skill score is equivalent to the model efficiency of Nash and Sutcliffe (1970). Defined this way, the skill score can be decomposed as the correlation coefficient squared minus two additional terms (Murphy 1988): where an r is the correlation coefficient between observations and predictions, σ is the standard deviation, and an overbar represents a time average. The second term, or conditional bias, is the difference in variance between the model and the observations and is zero when the slope of the regression line equals 1. The third term, or unconditional bias, is the mismatch of the means, and is equal to the intercept of the linear regression. The maximum skill is 1, and a skill of 0 represents a mean square error equal to the variance of the observations.

Combining Skill Assessments
Observational data used for model skill assessment can never have the spatial or temporal coverage of the model itself, and skills are only meaningful for the spatial and temporal coverage of the dataset used for the evaluation. If multiple sets of measurements are available, it can be useful to collapse multiple univariate skill assessments into single multivariate combined scores. Stow et al. (2009) presented the multivariate cost function, J, as a single measure of overall model skill: where x o and x p are the observations and predictions for all variables at all observed locations, T indicates the transpose of a vector, and R −1 is the inverse of the covariance matrix. The cost function simultaneously assesses model skill across multiple variables with disparate units, measurement frequency, and uncertainty. Note, however, that Eq. 3 provides equal weight to each observation. Combining multiple skill scores in this manner is most feasible when the observational data sets are similar, e.g., buoys that are deployed contemporaneously or CTD casts made at a common set of stations. When observational data are sparse, additional care must be taken in evaluating means across multiple dimensions to account for differences in the number of observations at different locations or times. O'Donnell et al. (2009) andMcCardell (2012) combined multiple Brier skill scores from distinct observational data by using mean relative errors which were then combined as a root-mean-square: where skill i are defined such that (1 − skill i ) represents the normalized model errors of the n individual skills to be combined. In Eq. 4, each individual skill (such as that from an individual dataset) receives equal weight in the overall score as opposed to each observation receiving equal weight. Any combination of skills must, of course, be judicious: individual skills by location, depth, season, or parameter contain far more information about model performance than any overall combined score.

Parameter Estimation and Improved Skill Through Data Assimilation
One method to increase skill is to use data assimilation to blend observations and dynamical models to make better estimates of parameters. The approach is generally valid when there is sufficient knowledge of the underlying processes' variability and sources of error. The differences between model and observations can be due to poor model representation of the physical or ecological process, deficient selection of model parameters, and/or observational errors. While data assimilation has often been applied to estuarine hydrodynamics (Spitz and Klinck 1998;Bertino et al. 2002;Xu et al. 2002;Frolov et al. 2009), the extension into ecological modeling is still mostly in development (Zhao et al. 2005;Cossarini et al. 2009). Data assimilation improves the statistical match between observations and the model but may result in physically unrealistic conditions elsewhere in the domain and can mask underlying deficiencies in the model. The Kálmán filter (Kalman 1960), which partitions observational and process uncertainties, is receiving increased application in ecology for its treatment of both time and space (Hinrichsen and Holmes 2009).

Considerations and Concerns Regarding Skill Assessment
No single approach to skill assessment is perfect, nor can one approach provide a comprehensive view of model performance. It is therefore important to use multiple lines of evidence in assessing model skill, for example, combining qualitative assessment of goodness of fit with quantitative skill metrics (e.g., Stow et al. 2009;Brush and Nixon 2010). Stow et al. (2009) also highlight that error exists in both observations and model predictions. Observational error results from sampling bias, analytical error, poor replication, and variability in both time and space. Prediction error results from imprecisely constrained parameters and numerical error. Given both sources of error, metrics that allow for both types of uncertainty will be most meaningful in assessing model skill.
Though traditional skill assessment almost exclusively focuses on comparison of observed and predicted states, recent work has highlighted the importance of accurately simulating key rate processes as well (Brush et al. 2002;Grangeré et al. 2009). This ensures that accurate model predictions of concentrations or biomass are achieved with reasonable inputs and outputs, i.e., that the model is Bgetting it right for the right reason.^For example, ecological models for prediction of hypoxia or other symptoms of eutrophication should accurately simulate not only phytoplankton biomass but also the rates of primary production and system respiration; this ensures accurate simulation not only of phytoplankton but also the production and consumption of oxygen, development of hypoxia/anoxia, cycling of nutrients, and flows of carbon. Coordination of observational and modeling campaigns in a given system is necessary such that process measurements and model simulations match in space and time.
Because skill metrics are typically based on model-data goodness-of-fit (Fitzpatrick 2009), which places emphasis on aligning the central tendency of a model with the central tendency of observations, Bimprovement^is often achieved through the addition of parameters that reduce model precision and increase model complexity. A skill comparison between a simple and complex model will almost always favor the addition of parameters because there is no penalty incurred from the increase in degrees of freedom. Burnham and Anderson (2002) use the Akaike information criterion to assess statistical model performance taking into consideration both model-data misfit and the number of parameters, but such an approach for complex mechanistic models is computationally impractical. Similarly, applications of numerical methods such as Markov chain Monte Carlo simulation and Bayesian techniques are emerging (Malve et al. 2007, Robson 2014). This approach is usually impractical in higherresolution coupled models, but Obenour et al. (2014) successfully applied Bayesian inference techniques to box modeling of hypoxia in the Gulf of Mexico. In any case, better knowledge about the expected statistical distribution of skill metrics could lead to optimization of the trade-off between accuracy and precision.

Balancing Spatial and Temporal Resolution
Hydrodynamic and ecological models are characterized by differing requirements for temporal and spatial resolution. Hydrodynamic models typically operate over time steps on the order of seconds with grid resolutions on the order of 10 1 -10 3 m. Ecological models may not need to operate at these same scales; typical biological rates like growth, respiration, and mortality fluctuate over periods of days to weeks, while other key rates (i.e., photosynthetic production) fluctuate over diel cycles in which hourly resolution is often sufficient. This is not to say that fine-scale biological interactions are unimportant to ecosystem-level processes, but current knowledge is insufficient to capture them reliably or to test hypotheses about their effects on system-level patterns. Thus, while fine-scale simulation may be necessary to adequately simulate ecological patchiness, the coarseness of the questions being asked and the limited availability of ecological data mean that ecological simulation at coarser scales may be appropriate for many model applications.
The trade-offs involved in the balance between spatial and temporal resolution in estuarine models is best illustrated through a review of model applications. In the following examples of model space-time complexity, Blow^spatial resolution is that greater than 1 km, and Blow^temporal resolution is a day or longer.

Low Spatial, Low Temporal Resolution
Despite recent growth in complexity and resolution, models of low spatiotemporal resolution continue to serve as valuable heuristic and management tools. A number of these are empirical models that describe the relationship between a variable of interest (e.g., chlorophyll, primary production, hypoxic volume) and one or more controlling variables (e.g., nutrient loading, chlorophyll; Cole and Cloern 1987;Nixon et al. 1996;Brush et al. 2002;Lee et al. 2013). Some of these empirical formulations (Evans and Scavia 2011;Forrest et al. 2011;Liu et al. 2011;Murphy et al. 2011;Turner et al. 2012;Lee et al. 2013;Scavia et al. 2013) are used operationally to produce annual forecasts of estuarine and coastal water quality in the Chesapeake Bay and northern Gulf of Mexico, with a focus on hypoxia/anoxia (IAN 2014;NOAA 2014). In general, these models predict variables at a single place or over a defined area and are primarily used to understand the dominant controls on a variable of interest and synthesize information across disparate study systems to look for general patterns. Such models do not include mechanisms explicitly, but mechanisms are often inferred based on our knowledge of the system derived from experimental and observational effort. Lower-resolution box models account for greater temporal and spatial variability than statistical models, but necessarily idealize transport and ecological processes. Nonetheless, these have successfully been used to analyze estuarine hypoxia (Hagy and Murrell 2007), quantify estuarine budgets (Boynton et al. 2008;Testa and Kemp 2008), and serve as ecosystem simulation models (Scavia et al. 2006;Ménesguen et al. 2007;Brush 2012;Lake and Brush 2015). Turner et al. (2006) developed a statistical model to understand and predict hypoxia in the northern Gulf of Mexico; this application represents a low temporal (annual) and spatial (O ∼ 10 4 km 2 ) resolution case study. The region has a welldocumented history of oxygen depletion in the coastal zone, ostensibly due to watershed inputs of nutrients from the Mississippi River. Extensive data of nutrient loadings (forcing data) and hypoxic zone area (prediction target) drove a statistical model that predicted over 80 % of the variation in hypoxic zone area over a 26-y period. The strongest drivers were nitrogen input and time; the time variable corresponded to temporally increasing carbon storage in benthic sediments, thereby increasing sediment oxygen demand. The authors used this mechanistic information to suggest that hypoxia could be mitigated by a focus on nitrogen load reduction rather than phosphorus. In this case, despite a lack of resolution in space, time, and process description, the statistical model proved its value in predicting hypoxic extent.

High Spatial, Low Temporal Resolution
Some ecological models focus on the description of ecological processes that do not change in a small period of time, but are spatially variable. An example is the presence or absence of seagrass which is temporally steady over a growing season but highly variable in space. Models with high spatial and low temporal resolution may address the influence of light availability and sea level rise on seagrass presence/absence (del Barrio et al. 2014), the effects of eutrophication, pollutants and macroalgae competition on seagrass biomass distribution (Giusti et al. 2010), seagrass shoot densities (Bearlin et al. 1999), and the conditions necessary for the restoration of submerged vegetation (Duarte et al. 2013). High spatial resolution is necessary for seagrass modeling as shoot density and root biomass vary with sediment type and light availability. Time resolution can be low (days, months, seasons, years) for seagrass dynamics due to low growth rates, long biological timescales, and seasonal variability of the different species (Alexandre et al 2008).
One example of a high-spatial/low-temporal resolution model is presented by del Barrio et al. (2014), where a coupled high spatial (10 m), low temporal (seasonally averaged) seagrass model was applied to West Falmouth Harbor, a eutrophic estuary in Massachusetts, USA. This model predicted seagrass presence/absence as a function of light availability with a spatial resolution of 10 m, whereas the temporal resolution was 2 months. This tool was used to evaluate future seagrass presence/absence based on nitrogen loading and sea level rise. With these two forcings applied separately, the model predicted an expansion of seagrass presence with nitrogen reduction whereas sea level rise reduced light availability to the canopy, reducing seagrass coverage. The combined effect suggested an overall expansion of seagrass meadows due to the reduction of nutrient loads, with sea level rise playing a secondary role over the short term. In the long term, once nutrient loadings reached a baseline level, sea level rise became more relevant as meadows migrated landward to adjust to the light climate. This offline, no-feedback model simplified the temporal evolution of seagrass meadows but resolved spatial variability in detail.

Low Spatial, High Temporal Resolution
Some model applications require that a process or system be simulated at high temporal resolution, but at low spatial resolution. The demand for this particular time/space trade-off results from the need to accurately simulate processes that operate over short timescales (e.g., organismal physiology), but without a substantial spatial representation. In estuarine environments, high spatial resolution is often required to capture hydrodynamic or bathymetric variability; if one is simply interested in algal physiology or plankton community interactions, a one-dimensional model may be sufficient, as with many applications of NPZD models (Fasham et al. 1990;Oguz et al. 2000;Soetaert and Middelburg 2009). Other examples of models with low spatial, but high temporal resolution include box models (Webster and Harris 2004) and sediment biogeochemical models (Hochard et al. 2010;Brady et al. 2013).
As an example of a low-spatial/high-temporal resolution model, Testa et al. (2013) applied a sediment flux model to simulate biogeochemical processes between the water column and sediment bed in Chesapeake Bay. The model spatially simplified the estuary to individual stations represented by three layers: the water column, the aerobic sediment layer, and the anaerobic sediment layer. The simulation time step was 1 h, which would be necessary to capture daily cycles of light, production, and respiration in applications where benthic algae were modeled explicitly. Application of this pointbased model indicated the importance of denitrification in the aerobic layer and illustrated spatial variability between multiple site locations. This type of model can be expanded horizontally (to a higher-resolution model domain) or vertically (with additional sediment layers), or it can be embedded in 3D coupled models to produce high-spatial, high-temporal resolution models.

High Spatial, High Temporal Resolution
From a realism perspective, simulation of estuarine hydrodynamics and ecology requires both high spatial and temporal resolution to capture physical and biological processes and their interactions. Many models that simultaneously simulate phytoplankton, dissolved oxygen, nutrient cycling, hydrodynamics, and sediment dynamics fall into this class (Xu and Hood 2006;Li et al. 2009;Zhang and Li 2010;Khangaonkar et al. 2012;Testa et al. 2014;Xue et al. 2014). Many models used in coastal water quality management opt for this highresolution configuration, and advances in computing power have removed most barriers to representing meter-scale bathymetry.
Examples of realistic, highly resolved models include 3D coupled hydrodynamic-biogeochemical models that have been used to understand the impacts of freshwater and nutrient loading on the biogeochemistry of coastal ecosystems (e.g., northern Gulf of Mexico, Fennel et al. 2011;Chesapeake Bay, Cerco and Noel 2013, Xu and Hood 2006, Testa et al. 2014Baltic Sea, Neumann and Schernewski 2008). Such models may simulate watershed and atmospheric loadings, hydrodynamics, sediment transport, water-column biogeochemistry, and sediment diagenesis and thus can be used to understand interactions between external forcing (river flow, wind, nutrient load) and internal biogeochemical cycling. With horizontal resolution of ∼1 km and time steps of <5 min, these models resolve bathymetric gradients and tidal timescale processes. In the case of the Chesapeake Bay and the northern Gulf of Mexico, model simulations have revealed the balance between stratification and nutrient loading in controlling anoxia (Cerco and Noel 2013), the spatial and temporal nature of the response of phytoplankton and hypoxia to elevated nutrient loading (Testa et al. 2014), and controls on phytoplankton biomass accumulation (Fennel et al. 2011). In all of these cases, the 3D coupled models provided understanding of the system beyond what previous statistical models provided, which illustrates the balance between generality (statistical model) and realism (fully coupled model).

Applications of Models Beyond Research
The models used for estuarine research and management span multiple spatiotemporal scales and levels of complexity. The size and complexity of the model as well as the software platform and user interface will reflect the purpose of the model and its intended user groups. In recent years, complex models that required a skilled modeler were favored by the management community. Examples include the USGS SPAR ROW model (Spatially Referenced Regression on Watershed Attributes; Moore et al. 2011), EPA WASP model (Water Quality Analysis Simulation Program; Ambrose et al. 1993), Chesapeake Bay Program Eutrophication Model (Cerco and Noel 2004), and Long Island Sound Study SWEM (Systemwide Eutrophication Model; HydroQual 1991).
While federal, state, and local managers still rely on complex models requiring scientific expertise, a growing trend has been to develop models accessible to the management community and other stakeholders without such expertise. In these cases, end users may alter model inputs within a reasonable range to assess the impact of various management decisions. A growing number of such models have been developed, especially in areas such as shellfish aquaculture where science, business, and management interact (e.g., Kellogg et al. 2014). The USGS has responded to the needs of the management community by providing a user interface for the SPARROW model that allows end users to explore the potential effects of changes in the watershed on the delivery of nutrients to streams and estuaries (Preston et al. 2011). Other models of nutrient delivery have been developed that allow users more control of locally relevant inputs, including the Nitrogen Loading Model (NLM, nload.mbl.edu; Valiela et al. 1997bValiela et al. , 2000 and the Regional Nutrient Management model (ReNuMa, Hong and Swaney 2013). Similar user-accessible models of estuarine response to nutrient loading and climate change are being made available to stakeholders through web browsers (Brush 2014). Improved user interfaces and a decrease in the technical skill required of the end user greatly expand access to a variety of model applications.
The advances in accessibility have also made it possible for models to be used in educational settings, from middle school classrooms to graduate level courses. While web-based, fastrunning models may be easily adapted for college courses, similar usage in secondary education requires cooperation between model developers and educators to translate Bcomplexm odels into easy-to-understand, interactive lessons. Ewing et al. (2003) noted that modeling has been Brelegated to the corners of the undergraduate curriculum^and that modeling has Benormous education potential but that potential is currently largely unmet.^Exploration of model results and participation in model building foster a more thorough understanding of ecosystem processes and naturally complement the analytical and quantitative skills emphasized in a STEMbased curriculum. The use of models for management decisions should be accompanied by a greater understanding in the general public of how models work, and integration of conceptual and mathematical modeling in secondary and undergraduate level programs presents that opportunity.
The utility of a model depends on a clear description of its attributes and limitations, including the assumptions and processes behind the model, potential errors and how error may be reduced, and the methods for verifying the results. Consideration of these issues is necessary for advancing the use of models in management and education, extending scientific knowledge into the realm of practical application. Potential obstacles can be navigated by (1) maintaining an open dialog with end users (via outreach workshops), (2) including relevant metadata in a standard format, (3) utilizing an interface that limits the users input to reasonable and appropriate simulations, and (4) providing accessibility and visualization tools that are intuitive.

Transparency and Utility: Making Model Output Accessible to End Users
Historically, numerical model code and output were opaque due to widely varying standards for documentation and longterm storage. Open-source philosophies and new web-based tools are shifting towards transparency. Modeling source code used to be modified in ad hoc fashion and then passed to other modelers via e-mail or FTP. Now, many source codes are stored in publically available revision control systems such as SVN or GIT. These systems allow users to browse code, easily track revisions, and recover the exact code used to perform simulations. This approach builds community understanding, an environment of trust, and the opportunity for reproducibility.
Model output used to be stored on investigators' machines and shared with users as derived products such as time series or section plots. Now, many models write self-describing machine-independent binary files such as NetCDF or HDF from which data can be efficiently accessed via standard web services, such as OPeNDAP (Cornillon et al. 2003). Software now exists to turn attributes stored in these files into ISO metadata records, allowing harvesting into data catalog systems designed to service a broad range of users such as data. gov. With data and metadata discoverable and accessible via services, scientists and developers can more effectively discover, access, and utilize model output (e.g., Signell and Snowden 2014). Data portals that leverage this infrastructure can then easily be constructed for targeted end-user communities. Online viewers, such as the MARACOOS asset viewer (assets.maracoos.org) and AOOS model explorer (http://data. aoos.org/maps/search/models-grids.php) enable non-expert users to visualize model output, map model results on top of different GIS layers, and compare with field data and with other models. In addition to serving output data sets, input datasets are also being shared via services, so users can examine input parameters and boundary conditions. Recent efforts to engage end users in participatory modeling have focused on encouraging a shared ownership and development of numerical tools that are designed to inform management. A recent example of this is BFishSmart,^a collaborative modeling exercise among academic modelers, commercial fishers, recreational anglers, and natural resource managers for the King Mackerel fishery (Miller et al. 2010). As has been suggested by social scientists (Paolisso et al. 2013), the success of such an effort rests heavily on the participatory process, in this case involving professional facilitators and a consensus-building approach to a series of workshops. Perhaps, the largest scale effort in this vein for restoration purposes has focused on the Everglades, where powerful economic forces have intersected with water quality goals and government management of water resources (van Eeten et al. 2002;Loucks 2006). In all of these instances, the modeling framework is the predominant decision-based tool being used to set regulatory actions (allowable catch, minimum catch size, nutrient loads, water use). Kelly et al. (2013) provide a decision tree to guide the selection of a modeling approach to be used for informing management, according to several considerations, including the reason for modeling, data availability, and the scale and processes of interest.

Recommendations
Modeling has emerged as a leading approach to the study and management of estuarine systems, and the importance of models will only increase during the current century. A variety of models exist that span multiple approaches (e.g., empirical, theoretical, mechanistic), spatiotemporal scales (e.g., seconds to years, meters to kilometers), degree of coupling (e.g., offline, online), and degree of complexity, and all provide important insights into estuarine system function. Models are working hypotheses about how a system functions, and no one model is Bright.^Similarly, no one model can singularly optimize the competing goals of generality, precision, and realism nor balance the trade-offs between complexity, over-parameterization, and uncertainty. In light of this review, we offer the following recommendations for future efforts in developing coupled hydrodynamic and ecological models of estuarine and coastal systems: (1) While the trend in modeling has been towards increasing resolution and complexity, more parsimonious, reduced, and intermediate complexity models offer useful alternatives that should be pursued in concert with complex, highly resolved models. These simpler models are particularly amenable to direct use by stakeholders and in STEM education. Modelers should consider the appropriate resolution and complexity necessary to address the research or management question at hand. (2) Modelers should rigorously evaluate model formulations and parameter values and test model output against theory and observations, rather than relying on Boff-the-shelf^models as truth. Increasing transparency in model code will allow the wider community (including model developers, empirical scientists, managers, and model end-users) to better evaluate model assumptions, formulations, and parameterizations.
(3) Training of the next generation of researchers should include cross-fertilization and development of skills in both observational and modeling techniques. Given the importance of developing physically and ecologically meaningful models rooted in empirical understanding, it will be increasingly important for modelers and empiricists to work together to develop more robust models. (4) Resources should be directed to develop and maintain long-term, spatially comprehensive observational databases, using accepted data standards. Given the importance of calibrating and independently verifying model performance against observations, these databases should be well documented and accessible. It is critical to verify model predictions of both state variables and key rate processes and to use a Bmultiple lines of evidence^approach that employs multiple assessments of skill. Research is needed to develop skill assessment metrics that allow for uncertainty in both model predictions and observations, allow for propagation of uncertainty in parameters, and account for the degree of parameterization. (5) Development of multi-model ensembles should be increased where possible. Application of multiple models, with varying spatiotemporal resolution and ecological/ biogeochemical complexity, may provide the best way forward in both heuristic and management applications of models. Comparison of models with varying physical and ecological resolution in the same system will provide important insights for future model development, and when applied to real-world problems, multiple models provide an envelope of likely response rather than relying on single, deterministic predictions. (6) Researchers should aim to translate state-of-the-art models into user-accessible, decision-support tools for managers and other stakeholders as appropriate. One promising way to do this is through participatory modeling approaches that integrate end-user feedback into model development. (7) Models should be integrated into educational activities to increase the public's fluency in modeling and serve as a recruiting tool for students into STEM fields. Models are ideal tools for demonstrating both scientific and mathematical principles, the critical linkage between these two fields, and the use of mathematics in real practice. Environmental models have the additional potential to put these lessons into a local context that simultaneously fosters environmental awareness and stewardship. Views, opinions, and/or findings contained in this report are those of the authors and should not be construed as an official U.S. Department of Defense position or decision unless so designated by other official documentation. Although the research described in this article has been funded in part by the U.S. Environmental Protection Agency, it has not been subjected to Agency review. Therefore, it does not necessarily reflect the views of the Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.