Toward a Comprehensive Water-Quality Modeling of Barnegat Bay: Development of ROMS to WASP Coupler

Located on the coast of New Jersey, the state with the highest population density in the United States (U.S. Census, 2010), the Barnegat Bay and Little Egg Harbor (hereafter Barnegat Bay) estuary is of important socioeconomic value to the region. The bay and its watershed attract many tourists, doubling the watershed population in the summer, and support a year-around commercial fishing industry with an overall $2 to $4 billion annual contribution to the New Jersey economy through tourism, employment, and ecosystems (Kauffman and Cruz-Ortiz, 2012). Increasing population and coastal development, however, pose a rising stress for this shallow back-barrier bay with limited exchange with the ocean.

For the last few decades, Barnegat Bay has experienced an increase in harmful algal blooms, oxygen depletion, and a degradation of habitat (Fertig et al., 2014; Hunchak-Kariouk and Nicholson, 2001; Kennish et al., 2007; Lathrop et al., 2001). The need for a comprehensive understanding of the relation between watershed land use and water quality in the bay required a bay-wide monitoring effort supported by an advanced modeling approach. For this purpose, the U.S. Geological Survey (USGS) New Jersey Water Science Center has partnered with the USGS Woods Hole Coastal and Marine Science Center to provide technical assistance and guidance to the New Jersey Department of Environmental Protection (NJDEP) in developing a total maximum daily load (TMDL) for the bay. Because of the complexity of the physics in the bay, coupling a three-dimensional hydrodynamic model with a water-quality model was proposed. The proposed modeling approach would also be advantageous for researchers and regulatory agencies in testing of various management alternatives to improve the water quality in the bay.

The Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005), an open-source ocean modeling system, has been successfully applied in many coastal, lacustrine, and estuarine applications (Banas, MacCready, and Hickey, 2009; Defne, Haas, and Fritz, 2011; Grifoll et al., 2013; Li et al., 2016; Matsumoto, Tokos, and Gregory, 2015; Ralston et al., 2013; White and Matsumoto, 2012). With the addition of the wetting and drying capability (Warner et al., 2013), ROMS is suitable for modeling tidal dynamics and water circulation, temperature, salinity, and sediment transport in a shallow, back-barrier estuary. The Water Quality Analysis Simulation Program (WASP; Ambrose et al., 1988; Di Toro, Fitzpatrick, and Thomann, 1983; Wool, Davie, and Rodriguez, 2003), developed by the U.S. Environmental Protection Agency (EPA), incorporates watershed nutrient loading and internal nutrient cycling to calculate dissolved oxygen (DO), oxygen demand, nutrient concentrations, sediment, and phytoplankton dynamics and has been widely used for water-quality and TMDL assessment in rivers, reservoirs, lakes, and estuaries (Abdelrhman, 2015; Camacho et al., 2014; Franceschini and Tsai, 2010; Hosseini, Chun, and Lindenschmidt, 2016; Kaufman, 2011; Lindenschmidt, 2006 and references therein; Tetra Tech, 2012, 2015). WASP can receive hydrodynamic information from other models (EFDC, DYNHYD, RIVMOD, CE-QUAL-RIV1, SWMM) through a binary hydrodynamic linkage file; however, a coupling between ROMS and WASP was not available prior to this study. Therefore, one of the main requirements of this study was developing a coupler between the two models, which can subsequently be applied to other studies and locations.

This paper presents the coupling methodology, together with preliminary results from the Toms River embayment, a subdomain of Barnegat Bay. Because of increased freshwater to saltwater transition and areas of higher flow, the Toms River embayment was selected to test the functionality of the ROMS-WASP Coupler (RWC) in a realistic application. First, the ROMS and WASP models and the RWC code are explained in the “Methods” section. This is followed by testing of the developed method at the Toms River embayment in the “Results” section. Potential improvements are discussed based on the results from this test case in the “Discussion” section. Finally, further recommendations and a summary are provided in the “Conclusion.”

Based on the prior experience and following earlier examples of successful coupling between WASP and other hydrodynamic models (Kim, Lim, and Cerco, 2011; Wool, Davie, and Rodriguez, 2003), a one-way offline (executed sequentially, where output from ROMS is passed to WASP for computation of variables) coupler was developed to study water quality of the bay. The hydrodynamic model setup was the same as in Defne and Ganju (2015) and spanned the entire bay (Figure 1a), whereas only the subdomain covering Toms River embayment (Figure 1b) was used to build the WASP grid (Figure 1c), which was built to demonstrate the functionality of the coupling with RWC. The details of the two models and the coupling between them are explained in this section.

ROMS is a free-surface, hydrostatic primitive equation ocean model with a terrain-following vertical and a boundary-fitted, orthogonal curvilinear horizontal coordinate system (Haidvogel et al., 2000). Finite differential equations are solved on a computational grid that is staggered in the vertical and horizontal (Figure 2; Arakawa C grid; Arakawa and Lamb, 1977), where vector quantities (e.g., horizontal and vertical velocity) are evaluated at the centers of cell faces (u, v, and w grid) and the scalar quantities (e.g., free-surface, density, and tracers) are evaluated at the centroid of cells (rho grid). Hence, one less node occurs in the u grid in x direction, and one less node occurs in the v-grid direction than number of nodes in the rho grid in each direction. Additionally, vertical-mixing variables are calculated at the bottom and top faces of each cell. This convention, while useful in some of the numerical computations in the model, requires caution while postprocessing the output, especially if switching between computer languages that treat arrays differently (e.g., from Fortran, where array indices start from 0 with [xi,eta] notation, to Matlab, where array indices start from 1 with [eta,xi] notation). For instance, for a computational rho grid with L+1 and M+1 total number of nodes in the x and y direction, respectively, the u-grid dimensions of [1:L, 0:M] in Fortran translates to [1:M+1, 1:L] in Matlab. The ROMS uses a Network Common Data Form (UCAR Community Programs, 2016) data structure for input/output.

The ROMS model setup, calibration, and skill assessment for the Barnegat Bay hydrodynamic simulation are described in detail in Defne and Ganju (2015). The same model was run for June–August 2012 to drive the WASP model. The hydrodynamic model was forced at the offshore open boundaries by tides (Mukai et al., 2002) and tidally averaged water level and current data from the ESPreSSO model (Wilkin and Hunter, 2013). Depth-varying salinity and temperature were also supplied at the open boundaries by the same model. Freshwater input was prescribed using USGS stream-gage data (USGS, 2012). Wind speed and direction, air temperature, and pressure, as well as longwave and shortwave radiation from the North American Mesoscale atmospheric model (NOAA/National Weather Service, 2012), were applied with a bulk flux parametrization as atmospheric forcing. The hydrodynamic model had seven evenly distributed vertical layers. The horizontal resolution varied between 40 m and 200 m in the bay. The Toms River embayment comprises 1569 cells per layer, with variable horizontal dimensions varying from 50 m to 160 m.

WASP is a dynamic compartment-modeling program for aquatic systems, including both the water column and the underlying benthos (Wool et al., 2003). The time-varying processes of advection, dispersion, point, and diffuse mass loading and boundary exchange are represented in the basic program. Application of WASP in estuarine settings involves specification of grid geometry and flow information, salinity, and temperature, which are obtained via linkage to a hydrodynamic model. Water-quality inputs to WASP are discussed below.

WASP was used to simulate water quality in the Toms River embayment during 1 July to 31 August 2012 after a 15-day model spin-up in June. This period was representative of eutrophication conditions in Barnegat Bay. In addition to the Toms River model, a full-bay scale water-quality model with fully simulated eutrophication processes is being developed separately. The full-bay model, which is intended for use in decision making, has gone through extensive calibration and validation, but its assessment has not been fully completed. Because the Toms River embayment presented in this study is used here for RWC demonstration purposes, providing any applied water-quality assertions or suggestions is beyond the scope of this study. Similarly, because of the extensive requirements in development and calibration of an advanced eutrophication model, a less complex eutrophication model was used in this study to simulate DO. The approach simulated DO, oxygen demands, nitrogen, and constant phytoplankton but did not include phosphorus and silica, particulate nutrients, inorganic solids, multiple phytoplankton groups, sediment diagenesis, settling, or advanced light extinction routines. Constant phytoplankton as chlorophyll-a (P-CHLA) concentrations that affected DO through photosynthesis and respiration were simulated, but the full range of diurnal variation in DO was not simulated.

The DO balance approach used is modified from Wool et al. (2003), as shown in Figure 3, where the arrows represent processes between constituents. For this study, simulated constituents included DO, carbonaceous biochemical oxygen demand (CBOD), ammonia (NH₄), nitrate (NO₃), dissolved organic nitrogen (DON), and P-CHLA. Particulate organic nitrogen (PON) was not simulated. Sources of DO include reaeration and phytoplankton photosynthesis. Sinks of DO include CBOD, sediment oxygen demand (SOD), and phytoplankton respiration. The CBOD is produced by breakdown of organic matter and is reduced to produce carbon dioxide (CO₂) through carbonaceous deoxygenation. Two classes of CBOD were simulated representing the watershed (CBOD 1) and bay water (CBOD 2) origin. Sources of DON include animal death/waste and algal death. The DON is converted to NH₄ through mineralization. NH₄ is oxidized to produce NO₃ through nitrification. The SOD is produced through sediment diagenesis and was specified as model input with this approach. Where observed data were not available, the calibrated parameters from the full-bay model based on the WASP advanced eutrophication module were used to implement the boundary conditions. Other inputs to the embayment model were derived from sampling data collected by NJDEP in watershed tributaries and the bay (Pang et al., 2017).

The WASP component of the Toms River embayment model comprises 38 segments per layer, with three total layers (Figure 1). Boundary loadings at the upstream end of the Toms River were derived from flow-weighted mean concentrations observed at in-tributary sampling stations BT03, BT04, and BT05 (shown in Figure 1, along with in-bay stations BB03, BB04a, and BB05a). Concentrations for north and south bay cross-sectional boundaries at the downstream end of the embayment were based on data collected at in-bay sampling stations BB03 and BB05a. Bay samples were collected at the surface and bottom for all constituents except NH₄, NO₃, and DON (calculated value), which were surface samples only. Initial concentrations in the model were approximated, as model spin up occurred in June. The BB04a was used to compare the model results with the measurements.

Calculated hydraulic and wind-induced reaeration values were simulated in the surface layer based on information from Covar (1976). Wind speed and air temperature were input for this calculation. The SOD rates, estimated from in situ benthic flux measurements by Wilson and DePaul (2016), were input in the bottom layer. Constant P-CHLA concentrations in each segment were derived based on average values during 16 June to 31 August 2012 from corresponding segments of the full-bay model. NH₄ and NO₃ loads were estimated for atmospheric deposition to the surface layer from data provided by the National Atmospheric Deposition Program (NADP; Schwede and Lear, 2014). NH₄, NO₃, and DON loads were estimated for unmonitored watershed drainages surrounding the embayment using simple regionalization methods based on hydrologic similarity (Sawicz et al., 2011). NH₄ and NO₃ loads associated with submarine groundwater discharge were estimated using a cell-by-cell water-budget analysis of existing groundwater flow model output (Cauller, Voronin, and Chepiga, 2016) and representative groundwater chemistry from near-bay discharge areas. Negative NO₃ loads are estimated to represent losses from west-side salt marshes of the bay because of denitrification.

Kinetic coefficients and rate constants used in the WASP model setup (Figure 3) are listed in Table 1. Limited calibration was done to improve the match between simulated and observed constituent concentrations by adjusting kinetics coefficients and constants. The main adjustments were to phytoplankton kinetics. The differences between ROMS and WASP models considered during development of the coupler and the details of coupling are explained next.

Because WASP simulates many eutrophication processes in each model segment and works on a single processor, the computational limit for the maximum number of segments is substantially smaller than that of ROMS, which uses parallel processing. Therefore, for large domains it is likely that a reduction in the number of grid members is necessary moving from a ROMS grid to a WASP grid. For this reason, the WASP grid was developed by aggregating the ROMS computational cells to larger WASP computational cells (segments). Aggregation presented congruent alignment of the WASP segment boundaries with the ROMS computational cells, thus minimizing the interpolations and providing more accurate mapping of volume and flow information between the two model grids than mapping between two independently structured grids.

The ROMS grid was transferred to a geospatial processing environment, where each computational cell is stored as a polygon feature. A much coarser draft WASP grid was overlaid with the ROMS grid to provide guidelines. The ROMS cells were aggregated to roughly an order of magnitude larger WASP segments by merging the cell polygons. A segment connectivity map displaying the neighboring information between WASP segments and a grid connectivity map displaying the index mapping between WASP segments and ROMS cells were created. Open boundaries, where nutrient loadings were to be assigned later during the simulations, were also defined at this step. A full scale Barnegat Bay model required mapping of 201,656 ROMS cells to 1827 WASP segments. For the Toms River test case, 4707 cells were mapped to 114 segments.

The RWC, which has been developed in Matlab for fast prototyping, is a one-way coupler (linker, or preprocessor) to aggregate ROMS fine output (in space and time) to the WASP coarser input. It provides mapping of ROMS computational cells to WASP computational segments, volume and flow information (volume, depth, velocity, and dispersion), and salinity and temperature. The suite of scripts developed for this purpose includes (1) segment information, (2) vertical mapping, (3) horizontal mapping, (4) time stepping, and (5) resampling modules (Figure 4). Each of these modules is briefly explained below and illustrated in Figure 5 with a transect from the Toms River embayment case. In the Toms River embayment case, each WASP layer has 38 segments (totaling 3 × 38 = 114 segments) and each ROMS layer has 1569 cells (totaling 7 × 1569 = 10,983 cells) in the horizontal, with an average of 96 ROMS cells per WASP segment. The exact number of cells aggregated per each WASP segment depends on the geometry of the segment in the horizontal and the number of layers aggregated in the vertical. In the vertical, the seven ROMS layers are aggregated to three in WASP. In time, the ROMS time step is 5 seconds (output every 30 minutes), and the WASP time step is 1 minute, retrieved by linear interpolation from ROMS output. The segment numbers in two adjacent layers are separated by the number of segments (38 segments). Numbering of flows start with the vertical flows, and top and bottom flow numbers (w1 and w2) are separated by the number of vertical flows per layer, which are also equal to the number of segments per layer. The horizontal velocity components start with an offset of total number of vertical flows (2 × 38 = 76 flows) and are separated by the number of horizontal flows (57 flows) per layer. Zero indicates a no-flow boundary (see Figure 1c and Figure 5).

While building the WASP grid, tables for describing interconnectivity of WASP segments and mapping of the ROMS cells to WASP segments are created with a GIS application. The first coupler module reads these two tables as input and stores the information in an array of SEGMENT data structure with a unique identifier (Segment.id) for each segment. For each segment, a field (center) to store the rho indices of ROMS cells contained in that segment, and fields (N, S, E, W) to store the list of rho-index pairs sharing a border at the interface with neighboring segments are created. Additionally, coordinates of segment centroid (centroid) to be used in calculating the distance between segments are computed. Segment boundary type (bc), indicating whether the segment has an open boundary, and the ROMS indices corresponding to that boundary (bc_eta, bc_xi) are assigned.

The second module implements the mapping between the two model grids and between the WASP segments in the vertical direction. An array of FLOW data structure is created in this module and is populated with the vertical components. The FLOW data structure comprises a unique identifier (Flow.id) for each interfacial flow; fields (Vxi, Veta) to store indices of the ROMS u-, v-, or w-velocity grids for flow through each WASP segment surface; a field (comp) that indicates flow component (u, v, or w); and a field (dir) that indicates the origin and destination segments of the flow. A field to store the distance (dist) between segment centroids in the vertical is created in this module but is populated later in the time-stepping module because the water-column height varies with time. The free surface at the top layer and the bottom surface of the bottom layer are set as no-flow boundaries. To relate SEGMENT to FLOW through each interface, SEGMENT structure is appended with vertical-flow fields w1 and w2 that store corresponding Flow.id values. Additionally, SEGMENT is also appended with indices of neighboring segments (topseg and botseg for top and bottom segments, respectively). Finally, the WASP layer of the segment (layerWASP) and the corresponding ROMS layers (layerROMS) fields are assigned. The vertical layer indices increment in opposite directions in WASP and ROMS. Therefore, if the top two layers of a seven-layer ROMS grid are mapped to the surface layer of a three-layer WASP grid, ROMS layers six and seven are mapped to WASP layer one (Figure 5).

In the third module, FLOW fields (Flow.id, Vxi, Veta, comp, dir, and dist) are appended with the horizontal flow information. First the interior, then the open-boundary horizontal flow components are assigned to these fields. This module appends horizontal flow fields (u1, u2 in W-E direction and v1, v2 in S-N direction) to SEGMENT. These fields, together with w1 and w2, point to the Flow.id of corresponding interfacial flows in FLOW and maintain integrity between the two structures. The SEGMENT is also appended with the neighboring segments information in all four horizontal directions (w, e, s, n for west, east, south, and north, respectively). Together with vertical neighbors, information from these fields is used in creating a segment map for reference. The segment map for Segment 11 is given as 11:w9,e14,s0,n90,b49,t0,Z1, indicating Segment.id, west, east, south, north, bottom, and top neighboring segments and the layerWASP, respectively.

The fourth module is the time-stepping module, where the output from the ROMS model is read at each time step and the corresponding fields of SEGMENT and FLOW are populated. The SEGMENT structure is appended with temperature (T), salinity (S), and segment velocity (segVel) fields, and the FLOW structure is appended with dispersion (E), flow rate (Q), interface area (A), and parameters crnu and brintt and iflowdir parameter fields. Segment velocity is used to calculate reaeration and/or volatilization rates but is not used to transport mass between segments. The segment velocity is calculated as the vector sum of interfacial flow-velocity components averaged at the centroid of the segment. The vertical dispersion coefficient is calculated as the average of vertical diffusivities of ROMS cells at the segment interface; the horizontal dispersion coefficient was set to zero. The coefficient crnu is the velocity between adjacent segments divided by the cross-sectional area between the segments and divided by the travel distance (distance between center of segments), and brintt is the dispersion coefficient times the cross-sectional area between the adjacent segments, divided by the travel distance. The field iflowdir is used to indicate the direction of flow (u, v, or w). Segment volumes and properties are calculated at the beginning of each WASP time step, and segment interfacial flows and properties are averaged during each WASP time step.

The fifth module is used for temporal resampling of the WASP input data. The maximum time-step size in WASP is limited by the minimum size of the WASP segments because of stability requirements. In addition to this, the frequency of the ROMS output interval itself is practically limited by the output file size. Therefore, resampling of the ROMS output at a higher frequency is required for linkage file. The SEGMENT and FLOW parameter values are resampled at a user-defined rate using three time steps centered about each time step (backward, center, and forward) from the original output. The interpolation scheme for resampling is selected by the user from various options (linear, cubic, or spline). Only the resampled time series between the center and forward time step are written out at each step. Finally, a Fortran executable (Wool, Davie, and Rodriguez, 2003) is called to convert this output to a binary hydrodynamic linkage file to be used as input for the WASP model.

The volume of each segment changes with time and is calculated by adding the volume of corresponding computational cells in ROMS. Flows through the segment faces are also calculated by aggregating the flows at the corresponding velocity grid points in ROMS. Because the horizontal dimensions of the WASP segments were usually an order of magnitude larger than the ROMS computational cells, it was also possible to run the water-quality model with a longer time step. Computational time steps ROMS and WASP were 5 and 60 seconds, respectively; however, because of the practical storage limitations, the hydrodynamic-model results were written to disk at a longer time interval, which required resampling of the intermediate output at a smaller time step as required by the water-quality model. In this case, ROMS output time step was 30 minutes but was resampled to 60 seconds with linear interpolation before inputting to WASP.

The RWC code can accommodate one-to-one or one-to-many connections (e.g., Segments 15 and 20 in Figure 1) between the WASP segments to address complex no-flow boundaries, large aspect ratios between adjacent model cells, or Courant stability criteria. During RWC tests, model results from one-to-one connectivity were compared to results from one-to-many connectivity for the same spatial area to ensure no additional mass balance errors were introduced.

Testing of the coupler was initiated using idealized one- and two-dimensional examples and comparisons with the EFDC model linkage with WASP. The idealized test cases were used as preliminary tests of correct transfer of information from ROMS to WASP and adequate connectivity between the two models. After these tests satisfied the flow and constituent mass balance requirements and produced rational results, RWC was tested with the Toms River embayment case. This test case is an estuarine subdomain with a complex, realistic, three-dimensional flow structure and therefore suitable to demonstrate application of the RWC. Evaluation of the RWC involved checking WASP input screens for correct passage of hydrodynamic information. The WASP Segments screen displays model geometry (i.e. spatial mapping between ROMS and WASP grids), segment volume, depth, and velocity. The WASP Exchanges screen displays vertical segment connectivity for dispersion. The WASP Parameters and Boundaries screens were checked to verify that they were populated correctly.

Water and constituent mass balance checks were important tests during the development of RWC. Water mass balance was a precursor to the constituent mass balance, hence the balance was checked by assessing volume-flow continuity in all model segments. The error introduced at each time step through coupling was quantified using the ratio of the volume change calculated from the volume of aggregated ROMS cells to the volume change calculated by the net flow through the interfaces of the WASP segment during that step:

where, Errⁱ is the percent absolute error at time step i, Vol is the calculated volume, Q_in and Q_out are calculated inflow to and outflow through segment interfaces, respectively. Mean and maximum errors for the entire simulation period were calculated with:

and

where, Err_avg is the mean error, Err_max is the maximum error, and N is total number of time steps. The mean error was minimized during the linkage development process until it was less than a targeted value of 1% in the majority of model segments. To evaluate constituent mass balance, a conservative tracer test was run in WASP using only the linkage file without any boundary concentration/loads or kinetics and bypassing all processes. Because WASP automatically sets boundary concentrations to a unit value (e.g., 1 mg/l) when importing a linkage file, concentrations should approach this unit value throughout the model domain if the simulation is run long enough.

Additional tests of the linkage file included comparing salinity and temperature output from WASP to that from ROMS. Advective flow output from WASP was evaluated to verify the tidal signature in the bay, and dispersive flow output was checked to verify constituent vertical stratification.

The RWC testing of the Toms River embayment was performed at all segments to ensure that the stability and mass and constituent conservation were satisfied through the entire computational grid; however, for brevity, Segment 11, which contains the measurement station BB04a, was used in this section to display test results. In general, the mass balance error was mainly attributed to the higher frequency resampling of the ROMS output but was well below the targeted value of 1% through the entire grid (Err_avg = 0.07% and Err_max = 0.67%). In Segment 11, the mean error was 0.03%, and the maximum error was 0.22%. This resulted in almost identical time series when volume time series output from WASP was compared to ROMS output calculated as the total volume of cells in Segment 11 (Figure 6a). Constituent mass balance was confirmed by introducing a conservative tracer at 1 mg/l concentration (default value in WASP) at the model open boundaries and by observing the concentrations asymptote to 1 mg/l in the entire domain after ∼15 days. This spin-up time was accounted for by starting the coupled models run 15 days before the start of targeted simulation period (i.e. 16 June). Measured surface salinity at BB04a was compared to the modeled salinity at the corresponding ROMS cell [k = 7, eta = 570, xi = 43] and the corresponding WASP Segment 11 (Figure 6b). Despite the spatial aggregation of 42 ROMS cells per layer and two ROMS layers in Segment 11, and resampling in time, the two models agreed well, with the difference between them being much smaller than the difference between the measured and hydrodynamically modeled salinity at BB04a. A transect along the Toms River embayment (Figure 1c) was selected to compare the salinity profiles. Two snapshots around the high- and low-tide times at the upstream boundary of Toms River are shown in Figure 7. Note that the ROMS cells along the transect were compared to the WASP segments created by aggregation (Segments 2, 4, 6, 7, 8, 10, 13, 16, 20, 23, and the segments below them). The salinity stratification was captured reasonably well despite the relatively much coarser WASP grid.

Simulation of trends in DO in the Toms River embayment is presented as a demonstration of the hydrodynamic linkage between ROMS and WASP. Simulated DO at site BB04a falls within the range of observed data (Figure 8) but exhibits less diurnal variation. Given that P-CHLA is applied as a constant value, phosphorus and silica are not simulated, and other eutrophication processes outlined previously are not included in the model; this result is not unexpected. Water temperature, which is passed through the linkage file, displays the expected relationship with DO: As water temperature increases, DO decreases, and vice versa. Salinity (not shown), which is also passed through the linkage file, exhibits a similar but less distinct inverse relationship with DO.

In terms of spatial distribution, simulated DO is higher upstream in the embayment than at the mouth (Figure 9). For example, the difference in average DO from surface segment 7 to surface segment 26 is 0.8 mg/L, and the difference between bottom segment 83 (below segment 7) and bottom segment 102 (below segment 26) is 0.48 mg/L. Based on 15 surface/bottom samples collected during July–August 2012, observed average surface DO is higher at BB04a than at BB03 and BB05a (7.19 vs. 6.68 or 6.30 mg/L, respectively). Observed average bottom DO is lower at BB04a than at BB03 and BB05a (5.79 vs. 6.46 or 6.13 mg/L, respectively). Additional calibration of reaeration, light extinction, and/or SOD could improve the agreement between simulated and observed DO in the lower layer.

In terms of vertical stratification, simulated DO is higher in the surface layer than the bottom layer both upstream and at the mouth of the embayment. For example, the difference between surface and bottom average DO in segments 7 and 83 is 0.48 mg/L and in segments 26 and 102 is 0.16 mg/L. Observed average DO is higher in the surface layer than in the bottom layer. Differences are 1.39, 0.22, and 0.16 mg/L at BB04a, BB03, and BB05a, respectively. Also, based on this information, simulated average DO exhibits greater vertical stratification upstream vs. at the mouth of the embayment. Observed average DO at BB04a likewise exhibits greater vertical stratification than at BB03 and BB05a.

Surface and bottom layers generally exhibit less DO deficit upstream than at the mouth of the embayment (Figure 9c,d). Surface-layer DO deficit occurs 7.8% of the time at segment 7 vs. 86.1% of the time at segment 26. Bottom-layer DO deficit occurs 36.2% of the time at segment 83 vs. 92.5% of the time at segment 102. Differences in DO deficit may be attributable to higher values of constant P-CHLA applied in the eastward direction, resulting in higher phytoplankton DO consumption. Alternatively, lower total nitrogen (TN) simulated at the mouth of the embayment may result in lower phytoplankton DO production.

The general trend of the observed data at site BB04a for each nitrogen constituent is reproduced by the model (Figure 10). The decline in simulated TN, DON, and NO₃ in late July/early August may be related to high flows from the Toms River to the embayment at that time. Because simulated nitrogen constituents match better with observed data than simulated DO with observed data, the inaccuracy of the DO match is likely less because of modeling of nitrogen processes. If calibration of DO was to be improved, emphasis would be on additional data collection to support the use of this approach.

During the development of RWC, a number of technical issues had to be resolved to minimize flow mass balance errors and to maintain conservative transport of the constituents. Moving from the hydrodynamic model to the water-quality model, reducing the spatial resolution of the computational grid has the potential drawback of loss of finer scale hydrodynamic information, and temporal upsampling introduces additional error, both of which affect the accuracy of constituent transport in WASP. In the final version, water mass balance error was reduced to less than the targeted value in all of the WASP segments, but an overall tendency toward larger error in shallower segments was noted. This issue requires further investigation with shallower bathymetry, especially for cases where frequent wetting and drying might occur in the domain.

The performance of the water-quality model appeared to be sensitive to the aggregation in the vertical direction. For example, the ROMS model had seven vertical layers that were collapsed to three layers in WASP by aggregating the top two layers in ROMS as the surface layer in WASP, which still resulted in a thin surface layer in some locations. Preliminary tests indicated that if the volume was decreased to a small value in these segments, constituent concentrations increase rapidly, causing model instability or run termination. Additionally, the calculated hydraulic reaeration rates in these small volume segments can be very high, making calibration of DO difficult (and insensitive) in the surface layer. Aggregating more ROMS cells into a WASP segment in the surface layer may help alleviate this issue.

The methodology and the scripts used in this study are applicable to other cases, provided the WASP grid is generated by aggregating ROMS computational cells. If future versions of WASP can accommodate more segments such that less aggregation of ROMS cells in WASP segments is needed, error related to reducing the spatial resolution may be reduced. Because of WASP stability requirements, however, smaller time steps are required, which would compound another issue: the size of the binary hydrodynamic linkage file, especially if the simulation period is also longer. For example, a 1-year long simulation for the full bay at a similar spatiotemporal resolution as the Toms River embayment case was estimated to require approximately 20 GB of storage space. Larger linkage files require longer runtimes in WASP, which prohibits numerous calibration attempts. Currently, no resolution has been suggested for this issue.

The development of the RWC to address the water-quality modeling needs of Barnegat Bay was inspired by previous successful examples of coupling WASP with hydrodynamic models and was motivated by its potential application to other geographical locations. The coupling required reducing the higher resolution of the hydrodynamic model grid and resampling of the hydrodynamic output at higher frequency. The developed code was tested for proper passage of hydrodynamic information between the two models and for satisfying the water and constituent mass balance. Constituent concentration results from WASP conservative tracer runs were evaluated to check for transport issues via absence of tidal signature and vertical stratification. A coupling approach that supported one-to-many connectivity between WASP segments was shown to work as well as a one-to-one connectivity and allowed for more flexibility in designing the WASP computational grid. The RWC provides the flexibility to use different spatial resolutions for hydrodynamic and water-quality models, which enables running only one higher resolution ROMS model to address physics-related questions of the system while coupling it to a WASP model at a lower resolution to answer water–quality-related questions. This can also facilitate building WASP models over existing ROMS applications. The ability to create and accommodate WASP grids with a one-to-many connection between the segments enables addressing any potential issues that may arise because of shallow depths or high Courant numbers by reshaping the problematic segments locally.

A WASP of the Toms River embayment was made for 1 July to 31 August 2012 to demonstrate the RWC. An intermediate-level eutrophication setup was used, which simulates only certain water-quality constituents, including DO, oxygen demands, N species, and constant P-CHLA. Model inputs were derived from water-quality data collected by NJDEP in watershed tributaries to Toms River and the bay. Where observed data were missing, the input is supplemented from a full-bay scale water-quality model that is being developed separately. The full-bay model has an advanced eutrophication module and has gone through extensive calibration and validation, but its assessment has not been fully completed. For this reason, the Toms River embayment has been used in this study as an example for a realistic application of the RWC. Providing any applied water-quality assertions or suggestions is beyond the scope of this study.

The calibration of the embayment model focused mainly on adjustments to phytoplankton kinetics. Simulated DO solubility was inversely related with the water temperature and salinity as expected. Spatial and vertical comparisons between simulated and observed DO were consistent at three sampling stations in the embayment. The simulated DO time series fell within the range of observed data at the nonboundary sampling station but exhibits less diurnal variation, which can be attributed to limitations of the setup. Simulated N concentrations match well with observed data at this sampling station. These results indicate that RWC performs adequately in an estuarine application. The full-bay model will provide additional verification of the coupling in a more diverse geographical setting (e.g., large tidal currents at inlets and through a more complete representation of eutrophication).

Funding was provided by the NJDEP and the Coastal and Marine Geology Program of the USGS. The author would like to thank Neil K. Ganju, who originally designed and supervised this study, and Robert Jr. Ambrose and Hugo N. Rodriguez for their input in early stages of developing the coupler code. The ROMS and WASP model outputs are available at USGS ScienceBase (USGS, 2017) and/or at http://geoport.whoi.edu/thredds/catalog/sand/usgs/Projects/BBLEH/run071tRX/catalog.html. For queries and questions about the RWC modules and source codes, contact the corresponding author, Z. Defne. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.