A Hydrodynamic and Sediment Transport Model for the Waipaoa Shelf, New Zealand: Sensitivity of Fluxes to Spatially-varying Erodibility and Model Nesting

Numerical models can complement observations in investigations of marine sediment transport and depositional processes. A coupled hydrodynamic and sediment transport model was implemented for the Waipaoa River continental shelf offshore of the North Island of New Zealand, to complement a 13-month field campaign that collected seabed and hydrodynamic measurements. This paper described the formulations used within the model, and analyzed the sensitivity of sediment flux estimates to model nesting and seabed erodibility. Calculations were based on the Regional Ocean Modeling System—Community Sediment Transport Modeling System (ROMS-CSTMS), a primitive equation model using a finite difference solution to the equations for momentum and water mass conservation, and transport of salinity, temperature, and multiple classes of suspended sediment. The three-dimensional model resolved the complex bathymetry, bottom boundary layer, and river plume that impact sediment dispersal on this shelf, and accounted for processes including fluvial input, winds, waves, tides, and sediment resuspension. Nesting within a larger-scale, lower resolution hydrodynamic model stabilized model behavior during river floods and allowed large-scale shelf currents to impact sediment dispersal. To better represent observations showing that sediment erodibility decreased away from the river mouth, the seabed erosion rate parameter was OPEN ACCESS

reduced with water depth. This allowed the model to account for the observed spatial pattern of erodibility, though the model held the critical shear stress for erosion constant. Although the model neglected consolidation and swelling processes, use of a spatially-varying erodibility parameter significantly increased export of fluvial sediment from Poverty Bay to deeper areas of the shelf.

Sediment Transport Models
Field experiments carry a high cost and are hampered by difficulties of observing water column sediment fluxes during energetic conditions such as floods and storms, except at discrete points served by deployed instruments. Numerical models based on the relevant processes for transport can be used to extrapolate point observations to continuous spatial scales, beyond the spatial and temporal coverage of field experiments. Here, we present a numerical model that complements a 13-month field campaign on the Waipaoa shelf, New Zealand.
Three-dimensional circulation and sediment transport models, such as the Community Sediment Transport Modeling System (CSTMS; [1]) resolve horizontal and vertical gradients, all of which can be important in the coastal ocean. The CSTMS has been implemented within the numerical hydrodynamic model ROMS (the Regional Ocean Modeling System; [2][3][4][5]). Although increased model complexity and resolution carries a heavier computational load, a three dimensional model was necessary to represent the complex bathymetry, bottom boundary layer processes, and river plume dynamics on the Waipaoa River continental shelf, New Zealand.
Many three-dimensional coastal sediment transport models have either neglected larger-scale currents or simplified them by using temporal and/or spatial averages to specify currents at the model's boundary, e.g., [6][7][8][9][10]. For example, a numerical model for Poverty Bay, the coastal portion of the Waipaoa Sedimentary System, which accounted for wind, wave, tidal, and river plume processes was developed by [6]. At the open boundaries, [6] accounted for tides and allowed disturbances to propagate through the boundary by using Chapman [11], Flather [12], radiation [13], and no-gradient boundary conditions for the free surface, two and three dimensional currents, and tracers, respectively. Recently, however, numerical models of continental shelf sediment transport have specified conditions along open boundaries based on estimates of coastal currents, temperature, and salinity from larger-scale, lower resolution models [14,15]. Like these examples, we build on previous efforts by nesting a finer-scale grid within a larger-scale hydrodynamic model, thereby accounting for larger-scale circulation patterns. For the event-driven Waipaoa shelf model, nesting not only allowed us to account for larger-scale currents, but was necessary to increase the stability of the model by reducing the reflection at the open boundary of sediment and freshwater from the river plume.
In many coastal environments, sediment fluxes are also affected by seabed erodibility, which can be defined as the amount of sediment available for entrainment into the water column at a given bed shear stress (see [16]). The treatment of erodibility is a distinguishing characteristic of cohesive and non-cohesive models (see [17,18]). For models of muddy cohesive seabeds, erosion typically depends on the seabed's critical shear stress, τ crit , and an erosion rate parameter, M, which regulates the rate of sediment resuspension (e.g., [1,17]; see section 2.4). Observations show that both parameters may vary with seabed porosity, the depositional history of the seabed, biological processes and other factors, e.g., [19][20][21]. For instance, recently-deposited sediments were easier to erode than material from consolidated, older seabeds in the York River estuary [16]. Based on seabed erodibility experiments, a bed consolidation scheme in which critical shear stress varied in time, space, and with depth into the seabed, depending on depositional history of the seabed, has been developed and implemented within numerical models [18,22]. Here, we developed a simpler parameterization that modified the erosion rate parameter to account for spatial variations in erodibility, based on seabed microcosm erodibility experiments (see section 3.6).

Study Site: Waipaoa River Continental Shelf, New Zealand
Located on an active tectonic margin and draining a small mountainous catchment, the Waipaoa River delivers material to the ocean primarily during floods [23,24]. The Waipaoa exports about 15 million tons of sediment annually, primarily through either gullying or landsliding, depending on riverine conditions [24]. This material is primarily mud, with a median grain diameter of 8.5 μm during flood conditions, and approximately 1% of the load is sandy bedload [23,25]. Because of the river's small catchment, rain storms induce flooding throughout the drainage basin, and delivery of sediment to the coastal zone typically coincides with energetic oceanic conditions [23].
Riverine sediments are delivered to Poverty Bay, an about 50 km 2 embayment that opens onto the continental shelf through a 10-km wide mouth (see Figure 1A). A counter-clockwise gyre driven by river discharge and the Coriolis force typically dominates currents in Poverty Bay [26][27][28]. On the shelf, water velocities during January 2010-February 2011 were primarily along-shore, but switched direction often, with an average current of 1.6 cm s −1 to the NE and a mean speed of 26.3 cm s −1 (data obtained from Hale, R. and Ogston, A., University of Washington (UW) [29]; tripod set-up described in [30]). Local winds, as well as larger-scale wind driven currents, southward travelling eddies, and coastally trapped waves likely drive water velocities [26,[31][32][33]. Surface wind and swell waves on the shelf have average periods of 9-10 seconds and significant wave heights of 0.8-0.9 m, although longer-period waves can reach the shelf from the Southern Ocean [28,34]. Wave-induced motion dominated bed shear stress calculations by an order of magnitude compared to current-induced stress [30,35]. During 2010 at a tripod deployed in 50 m of water near the Southern depocenter, bed stresses exceeded 0.15 Pa, a threshold for fine-sediment resuspension, for 46% of the deployment period (data from [30]). Over decadal and Holocene timescales, sediment accumulation on the shelf has occurred in two bathymetric lows to either side of Poverty Bay, but deposition is more variable over day-to month-long periods. Tripod observations and model estimates indicated that material is temporarily deposited in Poverty Bay following floods, and then, in the subsequent days to weeks, waves resuspend sediment and currents carry it to the shelf [6,28,35]. Observations of 7 Be activities also indicate that deposition over month-long timescales varies, depending on weather conditions [36,37]. For instance, 7 Be inventories from successive research cruises indicated that both erosion and deposition of terrestrial sediments occurred over different parts of the shelf during the January 2010 flood [37]. Over longer timescales, seismic profiles and 210 Pb radioisotopes (22.3 year half life) indicated two depocenters with maximum accumulation rates of ~1 cm year −1 occur located in bathymetric lows to either side of Poverty Bay, bordered by the coast and offshore anticlines ( Figure 1; [38,39]).

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.

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.

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 n u ,max cw  ,max cw  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][3][4][5]) and reduce the pressure-gradient truncation error [41][42][43][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.

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.

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].
,max cw 

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: (1) As indicated in Equation (1), the model assumed continuous deposition so that F cs,ised , the net entrainment of suspended sediment for each class, was calculated as the difference between the erosion of each sediment class, E ised , and estimated settling to the bed based on settling velocity of each sediment class, w s,ised , and on C s,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, f ised was the fraction of the seabed composed of sediment class ised, was the magnitude of the total wave-current-bed shear stress, and τ crit 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, z a , was specified as: where k 1 and k 2 were constants set equal to 0.007 m 2 s 2 kg −1 and 6.0, respectively, and D 50 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 (M max and M min ), and a transitional water depth (h transition ): ,1, , Based on observations, M max , M min , and h transition , were set to 4.5 × 10 −4 kg m −2 s −1 , 0.1 × 10 −4 kg m −2 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).

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][63][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: Nudging-Radiation OBC, evaluated at model boundary: (6) Here, ζ′ was the variable of interest (e.g., velocity or tracer concentrations) before nudging, T R,i and T R,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, 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.

Model Initialization and Forcing
Observed and modeled datasets used to initialize and force the model are listed in Table 3 and discussed in this section. All interpolations used linear Delaunay triangulation, unless specified below.

Model Grid Construction and Bathymetry
Designed to include the river mouth, Poverty Bay, and the proximal continental shelf, the model grid ( Figure 1) also encompassed the three depocenters identified by [23,38,39]. Use of a curvilinear horizontal grid and a stretched terrain-following vertical grid allowed the model to resolve regions of interest (e.g., two depocenters landward of shelf anticlines) and the near-bed and near-surface areas that have high vertical gradients in sediment concentration or velocity. While our model achieved lower horizontal resolution within Poverty Bay than that used by [6], it included more of the proximal continental shelf, allowing us to focus on shelf transport mechanisms.
The model grid had a horizontal resolution of about 450 m on the mid-shelf and was curved to reduce the number of terrestrial grid cells and to approximately parallel bathymetry ( Figure 1). This facilitated post-processing of data (i.e., across-shelf fluxes), and reduced model errors associated with nesting and along-isobath flow. In the unmasked (water) section of the grid, the angles of grid corners deviated from perpendicular by a maximum of 1.   [69]), while S. Kuehl (Virginia Institute of Marine Science; VIMS) provided multibeam data of the continental shelf and slope that had been obtained in 2005 on the R/V Kilo Moana [39]. S. Stephens (National Institute of Water and Atmosphere (NIWA), New Zealand) provided complete, though low resolution, 10 m bathymetric contours of the continental shelf and slope [70]. Finally, the historical gridded bathymetry from NIWA [71] provided the only data coverage near the entrance of Poverty Bay. The datasets were all referenced with respect to the WGS 84 datum and the universal transverse Mercator horizontal projection, except for the historical gridded dataset, for which the projection was unknown.  River observations [25]  Observed seabed properties [78]  ADV and OBS data [30]  Gust microcosm erodibility experimental data [36] Seabed characteristics for comparison to model estimates  Radiometric and X-ray analysis of cores [36,37] Comparisons revealed systematic offsets between the bathymetric datasets. In areas that overlapped (see Figure 3), the NIWA contours and Kuehl water depths were ~3 m and ~2 m shallower, respectively, than those from J. McNinch. The offsets were removed by adding 3 and 2 m to the NIWA contours and Kuehl data, which aligned them with the McNinch data from Poverty Bay. The deeper-water datasets were referenced to the McNinch data because bed stresses were most sensitive to bathymetry in shallow water and because 2-3 m was a smaller percentage of water column height in deeper areas compared to shallow areas. These three data sets were then combined for interpolation. Note that because the NIWA 10 m contours had lower resolution for much of the shelf, this dataset was relatively sparse compared to the gridded products. Thus, the gridded multibeam datasets dominated the model bathymetry where they provided coverage, while the NIWA contours filled in areas that were sparsely covered, primarily near the coastline and in the southwest and northeast portions of the grid.  Table 3). The National Institute of Water and Atmosphere (NIWA) contours also act as 10 m bathymetric contours. Note that the NIWA contours and multibeam data provided by S. Kuehl extend onto the continental shelf and slope.
After gridding, the model bathymetry was smoothed with a Shapiro (1975) filter to improve model stability [41,42,79]. Finally, water depths near the open boundaries of the model grid were adjusted to match the lower resolution ROMS-NZ bathymetry to facilitate model nesting. Both water depth and land-ocean masking were identical to ROMS-NZ over the region where nudging occurred. This methodology avoided noticeable seams where the datasets abutted, and the bathymetry, slope and curvature of the model grid was consistent with those of each individual dataset. The greatest uncertainty in bathymetry unfortunately lay at the entrance of Poverty Bay, where recent datasets on known projections did not provide much coverage.

Atmospheric Forcing and Waves
Modeled data were used as input to account for spatial and temporal variability in the wind and wave fields. Estimates from NZLAM and NZWAVE (described below; Figure 4c,e-f) were used for these because they had relatively fine resolution and were calibrated against local as opposed to global data. Values from these models, including wind velocities and wave estimates (significant wave height, mean wave direction, wave length, and mean wave period) were interpolated to the ROMS grid.
A local implementation of NOAA's Wave Watch 3 model [72], NZWAVE produced output every three hours on a 12 km resolution grid. Since NZWAVE did not provide wave spectra, the bottom wave period was assumed to be equal to the surface average period. This assumption resulted in underestimations of bottom wave period, because high frequency waves within a wave spectrum decay with water depth so that only the longer-period waves are felt at depth. Wave orbital velocities were calculated using linear wave theory and the interpolated NZWAVE data and bathymetry from the model grid. Linear wave theory was used instead of methods presented in [80] because estimates of wave orbital velocity based on this theory best matched tripod data.
Both model estimates and observed datasets provided atmospheric input. NZLAM, an implementation of the UK Meteorological (Met) Office's Unified Model [73], provided hourly estimates of wind velocity on a 12 km resolution grid. ROMS also required values for air pressure, cloud cover, precipitation, relative humidity, shortwave radiation, and air temperature as input. Hourly records of these meteorological data from the Gisborne, NZ airport were obtained from NIWA's National Climate Database web system [76] and applied uniformly across the Waipaoa shelf.

River Discharge
Waipaoa River water and sediment discharges were represented as a point-source entering Poverty Bay at a grid cell located at the river mouth. Observations of river stage were collected hourly at Kanakania Bridge, ~80 km upriver, above tidal influences, by G. Hall and D. Peacock at the Gisborne District Council (GDC) [77]. For water and sediment input, recently-calibrated rating curves provided by the GDC were applied to the river stage and water discharge, respectively (see Figure 4a,b). ROMS required that the vertical profile of the freshwater and sediment flux be specified at the point source. For the Waipaoa River, the profile was configured so that the freshwater and river sediment was delivered in the top half of the water column ( Figure 5).

Tides
Tidal velocities, amplitudes and phase components extracted from the Oregon State Tidal Prediction Software (OTPS) TPX07.1 global solution [74,75] were used to estimate tidal currents and sea surface height. OTPS accounted for eleven ocean tidal constituents, and was driven by satellite altimeter data (i.e., TOPEX/Poseidon and Jason). The sum of the OTPS estimates and values of water velocities and elevation from ROMS-NZ, which did not account for tides, were imposed at the model open boundaries.

Baroclinic Currents, Temperature and Salinity
Current velocities, temperature and salinity at and near open boundaries of the Waipaoa grid were nudged toward values from ROMS-NZ, a larger-scale baroclinic model adapted for northern New Zealand ( [65]; model used a similar model setup to [66]; see Figures 1 and 4d). ROMS-NZ provided estimates of current velocity, temperature and salinity every three hours on a two kilometer grid. This one-way nesting of the Waipaoa shelf model within the lower resolution hydrodynamic model allowed larger-scale circulation (i.e., shelf waves and offshore eddies) to influence modeled Waipaoa shelf hydrodynamics.
Three-dimensional, time dependent current velocities, temperature and salinity estimates from ROMS-NZ were linearly interpolated to the Waipaoa shelf grid and used for model initialization and nudging at model boundaries. ROMS-NZ estimates were unavailable for some grid cells near the coast in the interior of the grid where the land-ocean masking differed between the two models. At these sites, current velocities were initialized to zero, and initial temperature and salinity estimates were set equal to values from adjacent grid cells. Since land-ocean masking was identical between model grids near the open boundaries, these approximations only affected model initialization and not nudging near the open boundaries (see section 2.5).

Sediment Characteristics
Model calculations included a total of seven sediment types and eight seabed layers (Table 4). Different sediment classes were used to store fluvial and bed sediment so that model analysis could differentiate between materials from these two sources. Sediment classes were primarily distinguished based on their settling velocity, a primary control on transport in hydrodynamic-sediment transport models, e.g., [81], with nominal D 50 values provided based on observations from the river and seabed [25,78], and unpublished data [82] from the January 2010 cruise described in [36].
Estimates of effective settling velocity based on tripod measurements from the field site (obtained from A. Ogston and R. Hale, UW [29]; see [30]) informed our choice of settling velocities, while the initial seabed sediment distributions were based on sediment texture from Poverty Bay and shelf [36,69,78] (Figure 6). Observed grain size data were converted to percent sand and percent mud for all sites, and then spatially interpolated. Grain size data were unavailable for locations near the model open boundaries, so grid cells with water depth exceeding 300 m or near the northeastern and southwestern boundaries were assumed to be mud (sediment class 1). For the purpose of interpolation, grid cells at the coastline were assumed to be composed of sand (sediment class 2). Tripod-based estimates of effective settling velocities on the shelf ranged from <0.1 to ~1 mm s −1 , with a mode of ~0.1 mm s −1 during energetic shelf conditions at 40 m water depth, and were used to set the settling velocities of these two classes. Finally, the grain size distribution was constrained so that at least 10% of the seabed was composed of fast-settling material (sediment class 3) to enhance bed armoring, consistent with previous modeling efforts [6].  Four classes were used to represent sediment delivered fluvially. Their properties were informed by observations that estimated a median grain diameter of 8.5 μm in the Waipaoa River during floods [25], sediment budgets for the shelf, as well as the tripod estimates from the Waipaoa shelf (see above paragraph). Since information regarding the distribution of sediment settling velocity in the river plume was unavailable, a range of reasonable sediment settling velocities (0.15-1 mm s −1 ) was chosen. River mud was assumed to be flocculated because salinity profiles from research cruises in 2010-2011 [82] and results from other model runs [35] indicated that freshwater mixes quickly with ocean water, and very little stratification is observed in Poverty Bay. The fluvial sediment was logarithmically partitioned into these classes following sensitivity tests that considered the sediment budgets for Poverty Bay and the continental shelf following a three month model run representing early 2010. During this procedure, different settling velocities were prescribed for the slowest settling fluvial and seabed classes, 0.15 mm/s and 0.1 mm/s, respectively. This was not surprising, as the choice of fluvial settling velocities relied on an integrative view of sediment dispersal, i.e., the overall shelf sediment budget, while seabed settling velocities were selected based on local tripod observations on the mid-shelf ~10's of cm above the bed. Further study would be required to reconcile these, perhaps involving a model that includes particle aggregation and breakup.
Parameters related to erodibility (critical shear stress, erosion rate parameter) were informed by ADV and OBS (Optical Backscatter Sensor) measurements from the first two months of the tripod deployment (data from [30]) and Gust erosion chamber experiments conducted at the field site [36]. Critical shear stress for all sediment was set to 0.15 Pa, consistent with time series of tripod-based estimates of bed shear stress and sediment concentration, and Gust chamber erodibility experiments. For fast-settling grains (i.e., sediment classes 1 and 2), this critical shear stress was low compared to values based on the Shields parameter for the nominal grain diameters. However, in the absence of observations from which to specify various critical shear stress estimates for different classes of sediment, we used a uniform value. Note that calculated sediment fluxes were insensitive to the critical stress for fast-settling grains which settled to the bed quickly. Porosity was assumed constant across the shelf and set equal to 0.6 based on water content measurements [36]. Maximum and minimum erosion rate parameters in Equation (4) were assigned values of 4.5 × 10 −4 kg m −2 s −1 , and 1.0 × 10 −5 kg m −2 s −1 based on estimates from the Gust microcosm erodibility experiments [36]. Although the model neglected the role of depositional history and biology on seabed consolidation and erodibility, spatially-variable bed stresses and erosion rate parameters created a gradient of high to low seabed erodibility as discussed in section 4.1.

Sensitivity Tests
The Waipaoa shelf model described above, called the -standard model‖, was implemented for 15 January 2010-27 August 2010. Sections 4.2 and 4.3 evaluated the effect of model nesting and spatially-varying erodibility on model estimates using four sensitivity tests listed below. First, -moderately nudged‖ and -weakly nudged‖ test cases were run by increasing nudging timescales T R,I and T R0 from 2.5 days to 5 and 25 days to relax the degree to which larger-scale currents influenced model estimates. Second, -low M‖ and -high M‖ test cases applied spatially-uniform erosion rate parameters where M in Equation 4 was set equal to M min and M max instead of the spatially varying M used in the standard model. All sensitivity tests were run for the same seven month period as the standard model, except for the weakly-nudged simulation which became unsteady after about 20 March 2010. Sensitivity tests for other parameters were considered, but we chose to focus on this subset to demonstrate the importance of nesting for model stability, and because use of a spatially-varying erodibility has not been previously used with the CSTMS, yet provided a straightforward way to improve the agreement between modeled and observed patterns of erodibility and sediment retention.

Results and Discussion
Results for the standard model are evaluated, and then the sensitivity of estimates to model nesting and seabed erodibility parameterization is discussed.

Model Evaluation
Both the modeled and observed currents varied spatially and frequently reversed direction (Figure 7). For the most part, depth-and time-averaged current velocities on the shelf had a northeastward orientation, although the modeled and observed currents often had an offset in direction. At the 40 m deep tripod location, for example, velocities in the model were directed 15 degrees east of north, shoreward of measured velocities that were oriented 54 degrees east of north. Although the model underestimated peak water speeds at the three tripod locations, it replicated the spatial patterns of the tripod-observed time-and depth-averaged current speeds, which increased from the depocenter to the shallow site to the deep tripod ( Figure 8). Modeled velocities along the shelf-slope break were particularly fast, with the strongest shelf currents traveling along the shelf break and along isobaths, passing seaward of the southern anticline, through Poverty Gap, and inshore of the northern anticline. Lower current speeds were estimated in the lee of Mahia Peninsula, over the southern depocenter and in Poverty Bay. While shelf currents were generally to the northeast, a counterclockwise eddy formed within Poverty Bay, as seen in observations [21,27,28] and a previous model [6]. A persistent eddy also developed in the model over the southern depocenter, consistent with local tripod observations where depth-averaged currents fluctuated, but were primarily directed to the northwest or south. The modeled currents frequently reversed direction in response to boundary forcing from the larger-scale model. The reversals occurred on a timescale of days, similar to temporal behavior seen in measured currents (Figure 7). Discrepancies between the modeled and observed currents did occur, however. Factors that likely contributed to these include issues with the skill or resolution of the wind model; use of climatological temperature and salinity fields to force ROMS-NZ, which could affect the timing of large-scale eddies propagating southward along the shelf break; imperfect or smoothed bathymetric data; grid discretization; or choice of mixing schemes. Model estimates of waves, bed shear stresses, and sediment concentrations also captured the timing of observed episodic events (see [30]). Although wave periods were underestimated in the model, wave orbital velocities were comparable to values derived from tripod data (Figure 4f,g). Peaks in bed shear stress estimates occurred during the observed wave events, with a correlation coefficient of 0.72 and a bias of −0.20 to 0.02 Pa (Table 5; Figure 4h). Although time-averaged bed shear stresses were underestimated by about 17% for tripod deployments, the bias improved during periods of energetic waves and in the shallow area of Poverty Gap, where the largest sediment fluxes occurred. Modeled stresses were underestimated by about 5% at the 40 m deep tripod for time periods when ADV-derived data exceeded 0.15 Pa. For suspended sediment, we compared point measurements from a small sampling volume located 17-100 cmab, depending on tripod location and deployment, to estimates from the model's bottom grid cell, which represented a vertically-averaged value over the thickness of that layer (0.84 m at 50 m water depth). For these, the model underestimated near-bed concentrations from the tripod by more than 80%. Since the model cannot capture near-bed gradients in turbidity, the lower estimates for suspended sediment concentrations were not surprising. Consistent with observations, suspended sediment concentrations peaked during observed wave events (Figure 4i).  Erodibility in the model was evaluated by comparing estimates of seabed level variability to observations of eroded mass from Gust microcosm experiments (observations provided by [36]). Seabed level variability for each grid cell was calculated by taking the standard deviation of the time series of modeled seabed elevation (Figure 9). At each grid cell, seabed level variability depended on the erosion rate parameter, local sediment supply, and on active layer thickness (Equation (3); Figure 2) which varied with bed stress and grain size. Seabed level variability, as defined by the standard deviation, included both the net trend of erosion and deposition for a grid cell, as well as shorter-timescale fluctuations from individual floods and storms. For example, locations in Poverty Bay experienced relatively high rates of deposition that contributed to the high variability calculated there. Overall, seabed level variability was highest in areas shallower than ~30 m depth and where fluvial deposition occurred. The decrease in seabed level variability from Poverty Bay to the depocenters was consistent with estimates of eroded mass from Gust microcosm experiments (data provided by [36]), which indicated that the most erodible sediments were also found in Poverty Bay and the mouth of the bay (R 2 = 0.4; Figures 10 and 11). Differences between the model and observations could occur because Gust microcosm experiments were instantaneous measurements, while seabed variability was time-averaged, or because the model neglected processes such as bioturbation, and variations in parameters including critical shear stress (see section 4.3). Patterns of erosion and deposition estimated by the model have been evaluated using seabed observations of 7 Be inventories that indicate recent deposition of terrestrially derived material ( [36,37]; Figures 12 and 13). Overall, both the observations and model estimates exhibited high spatial and temporal variability [35][36][37]. Areas of deposition often occurred in close proximity to areas where little or no sedimentation was detected or estimated, and the shape of the footprint of recent terrestrial deposition changed between each research cruise. Both observations and model estimates generally showed enhanced deposition to either side of Poverty Gap, landward of shelf anticlines, although model estimates of deposition were located in shallower water than observed high 7 Be inventories. During every research cruise, radioisotope signatures were also high in some part of Poverty Gap, consistent with model estimates that showed sediment deposits there, particularly following periods of high discharge. Often, these high radioisotope signatures and model estimates of deposition were observed at ~40 m and ~30 m in Poverty Gap, respectively.   Overall, sediment fluxes were likely underestimated, but patterns of transport and deposition were consistent with observations. For instance, episodic and energetic waves dominated bed shear stress calculations and determined the timing of seabed resuspension, as seen in observations and other modeling studies (see above, section 1.2). Relatively slow water speeds and low sediment concentrations compared to observations would cause sediment fluxes to be underestimated, especially in the on-and off-shore directions. However, the model captured the spatial pattern of velocities, the frequent reversals of current direction, and the time-averaged currents; also, changes to open boundary conditions, as presented in this paper and during early model runs, affected the shape and location of the flood deposit for individual events, but had little effect over the entire 13-month model run, suggesting that time-averaged fluxes were consistent with observations. Finally, the skill of this model was comparable to others for the region. For the Poverty Bay model [6], the authors calculated biases and correlation coefficients of −5 to 0 cm s −1 and ~0.75 for wave orbital velocity, and −11 to −5 cm s −1 and −0.33 to 0.67 for current velocities, comparable to values in Table 5 for the shelf environment.

Model Nesting
Comparing behavior of the standard model to ones that used less rigid relaxation timescales showed that model nesting helped account for larger-scale currents, improved current velocity skill, and increased stability. Evaluations of model performance for this study considered water velocities because currents are important for sediment fluxes. A. Ogston and R. Hale, UW, provided time series of acoustic Doppler current profiler (ADCP) data from three locations ( [29]; see [30], Figures 1 and 7). Overall, stronger nudging to regional circulation did not significantly affect the mean current speeds (Table 5, Figure 7). At the 40 m tripod, time-and depth-averaged current speeds and associated standard deviations for 15 January 2010-20 March 2010 were 10.4 ± 6.3 cm s −1 , 10.3 ± 8.3 cm s −1 and 9.3 ± 7.7 cm s −1 in the standard, moderately-nudged and weakly-nudged cases, respectively. However, stronger nudging near the open boundaries increased the frequency with which along-shore currents switched direction, better matching observations (Figure 7). At the 50 m tripod in Poverty Gap, correlation coefficients between model estimates and observed depth-averaged along-shore velocities increased from 0.21 to 0.42 to 0.47 as the intensity of nesting increased from the weakly-nudged to moderately-nudged to the standard implementation (Table 5).
Model nesting also stabilized currents in areas near the open boundaries, reducing the reflection of the river plume at the grid's edge. Without nudging, the model failed within a couple of days because of excessively high water velocities (over 2 m s −1 ) at the boundaries near Mahia peninsula and the northeast coast as the river plume reflected off of the open boundaries creating a gyre within the domain. As expected, stronger nudging limited both the formation of the gyre and reflection at the open boundaries. To evaluate model behavior, the flux of freshwater through the open boundaries was estimated as: (8) where was the velocity perpendicular to an open boundary, S was salinity, and S 0 was the background salinity of 35.1 psu. During the January 2010 river flood (28 January-15 February 2010), the cumulative freshwater flux through the open boundaries increased from 0.12 × 10 6 m 3 to 0.23 × 10 6 m 3 to 1.12 × 10 6 m 3 when the nudging relaxation timescales, T R,i and T R0 , decreased from 25 days to 5 days to 2.5 days. These fluxes of water were equivalent to about 1%, 2%, and 11% of the freshwater discharge into the grid. Note that the nudging of tracers from cells near the open boundary but within the grid (Equation (5)) removed freshwater from the grid, removing more freshwater for shorter relaxation timescales. Therefore, the estimates of freshwater flux through the boundaries represent low estimates which removed approximately double the volume of freshwater from the standard model compared to the weakly-nudged simulation.
Mean current velocities were sensitive to the strength of nudging. During the first tripod deployment, for instance, mean currents in Poverty Gap at 40 m water depth changed from 2.1 to 3.4 to 4.1 cm s −1 , and the direction of mean water velocity changed from 104 to 72 to 54 degrees counter clockwise from east for the weakly-nudged, moderately-nudged and standard simulations. Current direction fluctuated frequently, however, so sediment dispersal remained relatively consistent among the different model runs, especially over timescales of months.
Therefore, the partitioning of sediment among different areas of the system (e.g., Poverty Bay vs. the rest of the shelf vs. off the proximal shelf) was relatively insensitive to the nesting scheme. Over the nine months of the model run, from 15 January to 25 August 2010, there was a 6% decrease in the sediment escaping the proximal shelf, a 6% decrease in sediment on the shelf, and a 13% increase in sediment in the bay for the weakly-nudged test case compared to the standard simulation. Similarly, a numerical modeling study [15] found that nesting increased sediment export from the Mekong delta front by <5%.

Sediment Erodibility
Choice of seabed erosion rate parameter (M in Equation (4)) influenced the amount of, and location of deposition on the shelf. In general, estimates of resuspension and sediment export from Poverty Bay to the shelf increased with the erosion rate parameter (Figure 12). Sediment fluxes in shallow areas were particularly sensitive to the choice of M due to increased sediment resuspension where bed stresses were high. Dispersal of slow settling material that remained suspended for relatively long times was also sensitive to M. In contrast, dispersal of sediment settling at 1.0 mm s −1 was relatively insensitive to M; differences in estimated sediment budgets were within 2% of each other for the cases of low, to spatially-variable, to high erosion rate parameters. For sediment settling at 0.15 mm s −1 , however, sediment export from Poverty Bay between 15 January and 27 August 2010 increased from 26% to 44% to 50% for the three cases. Despite the increased influx of sediment onto the shelf during the high M test case, however, less sediment was retained on the shelf because the high erosion rate parameter encouraged resuspension and resulted in more sediment export from the shelf compared to the standard model ( Figure 14). For the low, spatially-varying, and high M test cases, 7%, 14%, and 9% of sediment settling at 0.15 mm s −1 remained on the shelf, excluding Poverty Bay. Overall, spatially-varying erodibility increased deposition on the shelf relative to Poverty Bay, consistent with radioisotope-derived estimates of deposition on month long timescales [36]. Results from the standard model were most consistent with studies indicating that about a quarter of riverine material has remained on the shelf over decadal to century-long timescales [38]. Use of other seabed parameterizations for erodibility that account for bed consolidation and variations in critical shear stress, e.g., [18,22], could further increase sediment export from Poverty Bay following floods, and further strengthen the spatial trend of decreased seabed level variability in deeper areas of the shelf. For instance, some models account for the dependence of seabed erodibility of muds on depositional history such that the seabed's critical shear stress increase and decrease following erosional and depositional time periods, respectively [22]. Utilizing this type of seabed scheme would likely create areas of low critical stress in depositional areas following floods, such as Poverty Bay, enhancing sediment export in the days following high discharge events. However, this erodibility parameterization requires additional information about observed seabed critical stress profiles, is more computationally expensive, and has not yet been used for many realistically implemented sediment models (e.g., [83]).

Computational Concerns
Many decisions in the implementation of this three-dimensional numerical model required tradeoffs between desired accuracy and spatial resolution, and computational limits. The model had a total of 118 × 287 horizontal grid cells, each with 20 vertical water column layers and 8 vertical sediment bed layers. A total of nine tracer variables were included (salinity, temperature, and seven sediment classes), in addition to the momentum state variables. To provide estimates that overlapped with the Poverty Shelf field experiment, the modeled time period needed to span 13 months, from January 2010-February 2011, and provide estimates of state variables, including velocities, tracer concentrations, and sediment bed characteristics, every three hours for each grid point. ROMS has been parallelized using MPI (Message Passing Interface), which allowed us to run the model on VIMS' High Performance Computing (HPC) cluster using 48 nodes. The full 13-month model run required 9 days to run to completion. Some choices of model implementation significantly slowed the computations, including the MPData algorithm for horizontal advection of tracers, and the nudging of currents and tracers near the open boundaries. These components of the model were, however, important for model stability.

Summary
This project built on previous efforts by using a nested hydrodynamic-sediment transport model with spatially-variable erodibility to examine sediment fluxes on the Waipaoa Shelf. A three-dimensional sediment transport model accounting for a river plume, winds, waves, larger-scale currents, and tides was developed and implemented for the Waipaoa Shelf, New Zealand. These processes were represented using the ROMS-CSTMS framework in conjunction with locally-validated observed and modeled datasets described above. By varying horizontal and vertical resolution in the model, we focused on the area of interest and boundary layer processes while maintaining sufficient model efficiency. Sensitivity tests indicated that nesting helped to stabilize currents near the open boundaries, reducing the reflection of the river plume there, but variations in nudging did not notably affect sediment budgets for this implementation of the model. In contrast, a spatially-variable erosion rate parameter was needed to increase the export of material from Poverty Bay and retention of material on the shelf.