**3. 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.

#### *3.1. 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.4° and a mean of 0.13°. Vertical resolution varied with depth, and was about 0.40 m, 2.0 m, and 0.40 m near the surface, mid water-column, and seabed at 24 m water depth at the entrance to Poverty Bay. In deeper areas, resolution decreased so that, for example, surface, mid-water column and near-bed layers were 0.84 m, 5.3 m, and 0.84 m thick at 50 m depth.

Four datasets that each had a different focus provided the basis for the model's bathymetry. Multibeam was used to map Poverty Bay in 2005 and 2006 by J. McNinch (now at USACoE; see [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.


**Table 3.** Datasets used for model initialization and forcing.

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.

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.

**Figure 3.** Coverage of bathymetric datasets near Poverty Bay mouth. Datasets are labeled by source (see 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.

*3.2. 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.

**Figure 4.** Time series of observed and estimated weather conditions on Waipaoa Shelf. (**a**) River and (**b**) sediment discharge from Gisborne District Council, New Zealand; modeled wind speed from New Zealand Limited Area Model (NZLAM) (**c**) and current speed from ROMS-NZ (**d**) averaged over the model domain; (**e**) significant wave height; (**f**) bottom wave period; (**g**) wave orbital velocity; (**h**) total wave- and current-induced bed shear stress; and (**i**) suspended sediment concentrations. For panels (**e**)–(**i**), black lines indicate model estimates made for the grid cell nearest the tripod site (see Figure 1b), and grey lines indicate acoustic Doppler velocimeter (ADV) (wave height; bed shear stress) and optical backscatter sensor (OBS) (suspended sediment concentration) observations provided by Hale, R. and Ogston, A. (University of Washington; [30]) from the tripod at 40 m water depth.

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.

#### *3.3. 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).

### *3.4. 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.

**Figure 5.** Vertical profile of river input showing the partitioning of momentum, fresh water, and river sediment at the river mouth.

*3.5. 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).

#### *3.6. 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 D50 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<sup>í</sup><sup>1</sup> , with a mode of ~0.1 mm s<sup>í</sup><sup>1</sup> 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].


**Table 4.** Sediment classes and their characteristics

For all sediment classes: Critical Shear Stress: 0.15 Pa; Sediment Density: 2650 kg m<sup>í</sup><sup>3</sup> ; Porosity: 0.6.

**Figure 6.** Initial distribution of seabed sediment classes showing fraction of (**a**) mud and (**b**) sand.

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<sup>í</sup><sup>1</sup> ) 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<sup>í</sup><sup>4</sup> kg m<sup>í</sup><sup>2</sup> s<sup>í</sup><sup>1</sup> , and 1.0 × 10<sup>í</sup><sup>5</sup> kg m<sup>í</sup><sup>2</sup> s<sup>í</sup><sup>1</sup> 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.

#### *3.7. 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 *TR,I* and *TR0* 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 *Mmin* and *Mmax* 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.

#### **4. 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.

**Figure 7.** Time series of tidally- and depth-averaged water velocities from the Poverty Gap tripod in 40 m water depth. Observations (grey) and model estimates, including the standard model (thick red line), moderately-nudged (maroon), and weakly-nudged sensitivity tests (black).

#### *4.1. 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.

**Figure 8.** Map of estimated time-averaged depth-averaged current speed (shading; m s<sup>í</sup><sup>1</sup> ) and direction (black arrows). Long white arrows with blue outlines indicate observed current direction for tripod deployments. Black bathymetric contours indicate every 10 m.

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).


**Table 5.** Model evaluation statistics calculated (**A**) for the standard model for all tripod deployments with available observations; and (**B**) for depth-averaged currents from all three tripods' first deployments for all nesting sensitivity tests. This was the only time period when all models ran stably.

1 Correlation Coefficient; 2 Ratio of the standard deviation of the model estimates to that of the observations; 3 Difference between the mean of the model estimates and mean of the observations.

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*<sup>2</sup> = 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).

**Figure 9.** Modeled seabed level variability. Shading indicates log10 of the modeled seabed level variability, equal to the standard deviation of seabed thickness (cm) for each grid cell.

Patterns of erosion and deposition estimated by the model have been evaluated using seabed observations of <sup>7</sup> 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–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 <sup>7</sup> 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.

**Figure 10.** Observations of eroded mass from [36] for Gust microcosm erodibility experiments from January 2010 (**C1**), May 2010 (**C2**), September 2010 (**C3**), and February 2011 (**C4**). Grey bathymetric contours are every 10 m until 100 m water depth. Figure reproduced with permission from [36], copyright © Kiker, 2014.

**Figure 11.** Comparison of seabed level variability to observations of eroded mass from [36] for Gust microcosm erodibility experiments for specific sites. Multiple observations at a single site from were averaged.

**Figure 12.** Deposition per percent of riverine load following the March wave event for standard, moderately-nudged, low *M* and high *M* cases.

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<sup>í</sup><sup>1</sup> and ~0.75 for wave orbital velocity, and í11 to í5 cm s<sup>í</sup><sup>1</sup> and í0.33 to 0.67 for current velocities, comparable to values in Table 5 for the shelf environment.

**Figure 13.** Observed <sup>7</sup> Be inventories from (**A**) January, and (**B**) May 2010 research cruise. Figure reproduced with permission from [36], copyright © Kiker, 2014.
