Thank you for the review questions.

1. The reference to Holleman et al. (2013) alongside Chan et al. (2018) is made to point out the significance of numerical diffusion. The connection to Voronoi meshes is drawn by the second citation in the bracket, i.e. Chan et al. (2018). Voronoi meshes exhibit fewer acute polygon angles, spanned by cell vertices that are acutely remote from the cell centroid, which factor into numerical diffusion. Of course, flow aligned meshes are even better to tackle numerical diffusion but are here impeded by tidal dynamics bar the mesh itself becoming dynamic. Have changed "Furthermore, Voronoi approximations exhibit a reduction in numerical diffusion vs Delaunay meshes [4, 5]" to "In terms of numerical diffusion [4], Voronoi meshes exhibit a reduction compared to Delaunay meshes [5]." in the pending Arxiv revision.

2. Absent GPU acceleration, running this simulation along 5 other hydrodynamic simulations took one week on a Ryzen 9 7950X3D. That translates to about one day if all the CPU cores are dedicated to the simulation. Unboudedness over time is not a concern for numerical diffusion as it has a nivellating effect. In fact, its dampening effect can mitigate overshooting at the expense of accuracy.  

3. The algorithm has been validated with an MMS, method of manufactured solutions, which has been submitted separately for publication. The MMS was realized by oscillating the seabed to match the flow field to an analytical solution. The method can be used to validate the algorithm and mesh. Have included points mentioned in 3. and 4. below in the pending Arxiv revision.

4. The method is first order in space and time to attain high resolution meshes (Figure 6) to resolve waves while remaining efficient in terms of flops: To resolve waves, the cell size should be a log order below the part of the wave spectrum that is to be resolved. I.e. maximizing cell count and minimizing flops per cell are prioritized.

Timewise the LHS in equations such as #11 denotes forward Euler approximation. The ^{+\delta t} denotes a quantity at the subsequent time level. Past Delaunay versions of a species transport model (Lawen et al., 2013, Lawen et al., 2014 as cited in the paper), which worked in conjunction with other ocean models, offered for scalar quantities also semi-implicit matrix reordering algos. However, these attained only a tripling of time steps at the expense of flops for the reordering, rendering the net gain in terms of flops questionable. 

Based on this experience, a split approach is pursued: fast explicit Voronoi algorithm (small timesteps but less flops per step) and slow implicit reference solver (published separately) for algo cross-validation in addition to validation with MMS and survey. 

Yes, Voronoi schemes can be expanded to n dimensions. Interestingly, in this context, that might not improve results: for example, in the coastal case, the usual approach of resolving the vertical rather via multiple layers retains an alignment with the dominant horizontal current components and, thus, avoids numerical diffusion. That is, utilizing multiple layers achieves quasi flow alignment for the vertical. This example harkens back to question 1 above about numerical diffusion.

But perhaps expanding the Voronoi tessellation to the vertical could be an enhancement to model wave breaking and dynamic coastal meshes (4D Voronoi).

5. The model is capable of back-coupling modeled quantities onto hydrodynamic properties. There are two ways to  account for bed evolution: A.) Bed thickness change fluxes are modeled for representative periods, such as a neap-spring cycle, and then extrapolated for an annual seabed update. B.) To really simulate multiple years. The latter would preferably be supported by GPU acceleration.

6. Yes, did change "Numerical diffusion for flow-aligned unstructured grids" to "Numerical diffusion for flow-aligned unstructured grids with applications to estuarine modeling". Thank you.