*1.3. Objective*

Though both seabed and water column data have been collected for the Waipaoa River continental shelf, knowledge of sediment transport mechanisms is benefited by development of a three dimensional hydrodynamic-sediment transport numerical model providing spatial coverage unattained by observational efforts. This paper describes the implementation of the ROMS-CSTMS numerical model for the Waipaoa continental shelf and examines the sensitivity of sediment flux estimates to model nesting and seabed erodibility parameterizations.

**Figure 1.** Study site on North Island, New Zealand. (**A**) Waipaoa Shelf map showing locations of Poverty Bay, Poverty Gap, shelf depocenters (Dep), and the Lachlan and Ariel Anticlines (L.A. and A.A., respectively). Red arrow indicates river mouth. Grey bathymetric contours drawn every 10 m up to 70 m depth, while dashed grey line indicates shelf break at 150 m. Inset shows location of study site within the larger-scale model. (**B**) Waipaoa Shelf map showing shelf bathymetry up to 80 m water depth (shading; in meters), tripod locations (brown dots) and multi-core stations (black dots) from the January 2010 research cruise. Grey and black bathymetric contours drawn every 10 m to 100 m depth, then at 150 m depth. (**C**) Waipaoa model grid showing bathymetry, where each box shown encompasses 25 model cells. Model boundaries labeled "NW", "SW", "SE", "NE". (**D**) Bathymetric map of larger-scale model Regional Ocean Modeling System (ROMS)-NZ covering the eastern half of the North Island, New Zealand. Black box indicates location of Waipaoa shelf model.

#### **2. Model Development**

This section describes the equations and numerical schemes used to specify hydrodynamic and sediment transport processes within the model and at the boundaries of the grid. Table 1 lists symbols for all equations.

#### *2.1. Hydrodynamic Model and Numerical Schemes*

ROMS-CSTMS, a community-developed numerical circulation and sediment transport model, solves the equations for Reynolds-averaged Navier-Stokes, tracer advection-diffusion, and continuity using the hydrostatic and Boussinesq assumptions as described in [3]. Tracer concentrations could represent an array of different tracers, which included temperature, salinity, and seven sediment classes. River discharge was treated as a point source for momentum, temperature, salinity and suspended sediment. Sources and sinks within the governing equations included, but were not limited to, bottom friction, wind stress, and nudging to match the regional grid (see section 2.5). The density equation of state accounted for temperature, salinity, and sediment concentrations [1,40].

ROMS distinguishes itself from other community hydrodynamic models by its model grid, and time-stepping and advective schemes. It uses a curvilinear orthogonal grid in the horizontal and a stretched, terrain-following grid in the vertical which allows it to carry high resolution in both the surface and bottom boundary layers [2]. The numerical schemes in ROMS include split barotropic-baroclinic modes (Leap-Frog—Adams-Moulton predictor-corrector scheme; [2–5]) and reduce the pressure-gradient truncation error [41–44] by redefining the pressure-gradient term [3,45,46]. ROMS also provides high-order schemes for estimating both vertical and horizontal advection [3,47]. For advection of sediment and other tracers, ROMS provides the MP-DATA scheme (Multidimensional Positive Definite Advection Transport Algorithm; [48]) to avoid numerical oscillations and negative concentrations, and reduce numerical dispersion, e.g., [1,49]. For vertical advection of sediment, the ROMS framework implemented the PPM (Piecewise Parabolic Method) so that relatively large timesteps can be used for faster settling sands without the introduction of instabilities [1,50,51]. These numerical schemes and reasonably high spatial resolution are important for representing the high gradients typical of coastal and estuarine settings without sacrificing computational efficiency, e.g., [49]. For application of ROMS to the Waipaoa shelf, numerical schemes and the model grid (Figure 1) were chosen to reduce error in numerical computations without sacrificing model efficiency. A timestep of 15 s, and the numerical schemes listed in Table 2 were used.

#### *2.2. Surface Boundary Formulation*

The surface boundary formulation in ROMS was adopted from the physically-based COARE (Coupled Ocean Atmosphere Response Experiment) framework [52,53]. In this one layer boundary model, wind and rain transfer momentum from the atmosphere to the ocean.


**Table 1.** Model parameters as described in this paper.


**Table 2.** Numerical Schemes for Waipaoa Shelf model.

#### *2.3. Bottom Boundary Layer Formulation*

This implementation of ROMS used the Sherwood, Signell and Warner [1] bottom boundary layer parameterization, a physics-based approach that could account for form drag and ripples. This formulation, defined as "SSW" within ROMS, was based on [54] that divided the bottom boundary layer into two sections: a thin combined wave-current boundary layer, and a current boundary layer. For this study, the hydraulic roughness represented the grain size roughness, which was set equal to 0.30 mm based on estimates of bed shear stress from acoustic Doppler velocimeters (ADVs) provided by [30]. The model used the bottom roughness, the eddy viscosity profiles, wave orbital velocities provided by input files, and currents 20 cm above the bed, estimated by the hydrodynamic model, to calculate the total maximum current-wave-induced bed shear stresses, , following [54]. W *cw*,max r

#### *2.4. Seabed Model*

As summarized in section 2.1, the model accounted for suspended transport, erosion and deposition. Both fluvial discharge and seabed erosion provided sediment to the water column. To calculate erosion and deposition, Equations (1) and (2) were calculated for multiple sediment types, each having assigned values for settling velocity and diameter. Other values (e.g., erosion rate parameter, critical shear stress, and seabed porosity) were identical for all sediment classes. Processes not explicitly represented in the model include flocculation [55], seabed consolidation [56], and bioturbation [57].

Erosion was calculated following the Ariathurai and Arulanandan formulation [58] for each sediment class with index *ised*:

$$F\_{cs,load} = E\_{load} - \frac{\partial \mathcal{W}\_{s,load} C\_{s,1,load}}{\partial \mathbf{z}\_s} \tag{1}$$

$$E\_{load} = \begin{cases} M(\mathbf{l} - p) f\_{load} \left| \frac{\mathbf{r}\_{cve,\max}}{\mathbf{r}\_{cve}} \right| \cdot \tau\_{cvi} & \text{when } \left| \frac{\mathbf{r}}{\mathbf{r}\_{cve,\max}} \right| > \tau\_{cvi} \\ \mathbf{0} & \text{when } \left| \frac{\mathbf{r}}{\mathbf{r}\_{cve,\max}} \right| < \tau\_{cvi} \end{cases} \tag{2}$$

As indicated in Equation (1), the model assumed continuous deposition so that *Fcs,ised*, the net entrainment of suspended sediment for each class, was calculated as the difference between the erosion of each sediment class, *Eised*, and estimated settling to the bed based on settling velocity of each sediment class, *ws,ised*, and on *Cs,1,ised*, the mass of suspended sediment per unit area of the bed in the bottommost grid cell for each sediment class. In Equation (2), *p* was porosity, or void fraction of the seabed, *fised* was the fraction of the seabed composed of sediment class *ised*, was the magnitude of the total wave-current-bed shear stress, and *IJcrit* was a constant critical bed shear stress. Erosion during any time step was limited to the amount of each size class available within a thin active layer whose thickness, *za*, was specified as: *cw*,max W r

$$\tau\_a = \max\left[k\_1(\left|\tau\_{c\circ,\max}\right|-\tau\_{c\circ i}), 0\right] + k\_2 D\_{\circ 0} \tag{3}$$

where *k1* and *k2* were constants set equal to 0.007 m2 s2 kg<sup>í</sup><sup>1</sup> and 6.0, respectively, and *D50* was the median grain size on the seabed [59]. As implemented for the Waipaoa shelf, active layer thicknesses on the mid-shelf rarely exceeded ~5–10 mm and increased in shallow areas and near the anticlines due to the relatively high bed stress and the larger sediment grains found there (Figure 2).

Consistent with observations of erodibility on the Waipaoa shelf, the model formulation was modified to encourage erosion of sediment from shallow areas by varying the erosion rate parameter, *M*, with water depth, *h*. Choice of *M* used for a given water depth depended on parameters including the maximum and minimum erosion rate parameters (*Mmax* and *Mmin*), and a transitional water depth (*htransition*):

$$M = \left| \left[ \left( M\_{\text{min}} - M\_{\text{max}} \right) \Big| \left( 1 + \exp \left( 0.2 \bullet \left( h - h\_{\text{reaction}} \right) \right) \right) \right] + M\_{\text{min}} \right| \tag{4}$$

Based on observations, *Mmax*, *Mmin*, and *htransition*, were set to 4.5 × 10<sup>í</sup><sup>4</sup> kg m<sup>í</sup><sup>2</sup> s<sup>í</sup><sup>1</sup> , 0.1 × 10<sup>í</sup><sup>4</sup> kg m<sup>í</sup><sup>2</sup> s í1 , and 30 m, respectively (see section 3.6; [36]). Sections 4.1 and 4.3 evaluate the erodibility parameterization by comparing model results that used Equation (4) with two cases that used spatially-uniform *M* (see section 3.7). Both the active layer thickness and spatial variation in the erosion rate parameter influenced seabed erodibility. Other parameters such as critical shear stress likely also varied spatially and affected erodibility. However, estimations of parameters from erodibility experiments and tripod data can carry substantial uncertainty, and so this paper focused on a single variable, the erosion rate parameter, because it required no additional information about the seabed critical stress profiles and was computationally efficient. Use of a spatially-varying critical shear stress would likely also encourage erosion of sediment from shallow areas, similar to the parameterization used here. Alternate parameterizations from the literature are discussed in section 4.3.

Sediment bed properties such as grain size distribution were stored for eight seabed layers that each initially represented 20 cm of sediment. Erosion and deposition of multiple sediment classes modified the thickness of seabed layers and the grain size distributions stored for the sediment bed, as described in [1]. These changes impacted the upper few layers of the sediment bed, while the deeper layers served as a repository of sediment. The surficial seabed layer was ~1 cm across the shelf on average, but was thinner in areas of low deposition and where the active layer was thin (Figure 2).

**Figure 2.** Shows (**a**) time-averaged active layer thickness in the standard model, and (**b**) the spatially-varying erosion rate parameter, *M*, from Equation (4).

#### *2.5. Open Boundary Conditions*

The Waipaoa shelf model grid was bounded by land on the northwestern side (Figure 1), so a free-slip wall condition was used there which specified a zero gradient condition for tracers and sea surface elevation, set water velocities normal to land equal to zero, and used a free-slip condition for tangential velocities. Along the other three edges, open boundary conditions (OBCs) for sea surface height, barotropic and baroclinic velocities, and tracer concentrations accounted for tides, shelf waves, and the transient behavior of the river plume. Radiation conditions along the southwest, southeast, and northeast boundaries allowed waves to propagate through them without reflecting. Specifically, the Chapman [11] and Flather [60,61] conditions were applied there for the free surface and barotropic momentum, respectively, to account for tides [62–64]. Velocity and sea surface height at the boundary were required as input and were specified using data from ROMS-NZ, a larger-scale hydrodynamic model ([65]; model used a similar setup to [66]; see Figure 1; section 3.5). Sediment concentrations, which were not provided by ROMS-NZ, were nudged toward zero at the boundaries. This parameterization implied that external sediment sources were negligible, a common assumption in sediment budget calculations, e.g., [38,39], and in many models of riverine-dominated systems, e.g., [6,15]. Nudging sediment concentrations toward zero also assumed that material leaving the grid did not reenter the model domain, which was a reasonable expectation because the largest fluxes occurred during floods, when the river plume carried sediment off of the proximal

shelf [30,35]. Using the oblique radiation condition for baroclinic velocities and tracer concentrations reduced artificial reflections at the boundaries [11,13,67].

Similar to other studies [13,68], we nudged baroclinic current velocities, salinity, temperature, and suspended sediment concentrations within 30 grid cells of the open boundaries toward values specified from ROMS-NZ, or zero, as described above:

Nudging, evaluated within grid interior:

,

*R b*

*T*

*RO*

*RO OBC*

$$
\xi = \xi' + \frac{F\_{gvd}}{T\_{R,l}} \left(\xi\_{ohc} - \xi'\right) \Delta t \tag{5}
$$

Nudging-Radiation OBC, evaluated at model boundary:

$$\frac{\partial \mathcal{L}}{\partial t} + c\_x \frac{\partial \mathcal{L}}{\partial x} + c\_y \frac{\partial \mathcal{L}}{\partial y} = \left(\frac{1}{\mathcal{T}\_{\mathbf{x}, \phi}}\right) \tilde{\varphi}\_{\text{ac}}(t) \tag{6}$$

$$c\_x = -\frac{\partial \mathcal{L}}{\partial t} \frac{\partial \tilde{\varphi}\_{\text{ax}}}{\partial \tilde{\varphi}\_{\text{Ox}}^2 \dot{\lambda}^2 + (\tilde{\mathcal{C}} \frac{\partial \mathcal{L}}{\partial \mathbf{y}})^2}$$

$$c\_y = -\frac{\partial \tilde{\varphi}}{\partial t} \frac{\partial \tilde{\varphi}\_{\text{Ox}}}{\partial \tilde{\varphi}\_{\text{Ox}}^2 \dot{\lambda}^2 + (\tilde{\mathcal{C}} \frac{\partial \mathcal{L}}{\partial \mathbf{y}})^2}$$

$$= \begin{cases} T\_{RO} & \text{current velocities directed} & \text{out of} & \text{grid} \\ & T\_{RO} \text{ current velocities dominated} & \text{current width} \end{cases}$$

Here, *ȗƍ* was the variable of interest (e.g., velocity or tracer concentrations) before nudging, *TR,i* and *TR,b* were the relaxation timescales in the interior of the grid and at the open boundaries, and *ǻt* was the timestep. The larger-scale model, ROMS-NZ provided *ȗobc(t)*. The open boundary relaxation timescale, *TR,b*, decreased with decreasing *TRO* and increasing *FOBC*, enhancing nudging when currents flowed into the model grid. Based on sensitivity tests, *TR,i*, *TRO* and *FOBC* were set to 2.5 days, 2.5 days and 2.5. Section 4.2 presents results for "moderately-nudged" and "weakly-nudged" sensitivity tests where *TR,i* and *TR0* were doubled and multiplied by a factor of 10, respectively, increasing the influence of larger-scale currents. The parameter *Fgrid* specified spatial variability for nudging so that it was non-zero only for locations within 30 grid cells of open boundaries, and weakened sinusoidally with distance from the open boundary:

*T F current velocities directed in to grid*

° ° ¯ ¿

$$F\_{grid} = \begin{cases} 1 - \sin\left(\frac{\pi}{2}\frac{j}{30}\right) & for \quad j \in [1, 30] \\ 1 - \sin\left(\frac{\pi}{2}\frac{J - j + 1}{30}\right) & for \quad j \in [J - 29, J] \\ 1 - \sin\left(\frac{\pi}{2}\frac{I - i + 1}{30}\right) & for \quad i \in [I - 29, I] \\ 0 & otherwise \end{cases} \tag{7}$$

where *I* and *J* were the total number of grid cells in the NW-SE and SW-NE directions and coordinates *(i, j)* indicate location within the grid. Note that since ROMS-NZ did not include sediment, this meant that suspended sediment concentrations within 30 grid cells of the boundaries were nudged toward zero. Nesting not only enabled the model to account for larger-scale currents, but also reduced reflections of river plume salinity and suspended sediment concentrations at the boundary, which increased model stability.
