Coupled Delft3D-Object Model to Predict Mobility and Burial of Munition on Sandy Floor

Coupled Delft3D-object model has been developed to predict object’s mobility and burial on sandy seafloor. The Delft3D model is used to predict seabed environment such as currents, waves (peak period, significant wave height, wave direction), water level, sediment transport, and seabed change, which are taken as the forcing term to the object model consisting of three components: (a) object‘s physical parameters such as diameter, length, mass, and rolling moment, (b) dynamics of rolling cylinder around its major axis, and (c) empirical sediment scour model with re-exposure parameterization. The model is compared with the observational data collected from a field experiment from 21 April to 23 May 2013 off the coast of Panama City, Florida funded by the Department of Defense Strategic Environmental Research and Development Program. The experimental data contain both objects’ mobility using sector scanning and pencil beam sonars and simultaneous environmental time series data of the boundary layer hydrodynamics and sediment transport conditions. Comparison between modeled and observed data clearly show the model capability.

. Flow chart of the coupled Dlft3D-object model to predict objects' mobility and burial.
At the same time, a range of surrogate munitions were deployed that include variations in caliber, bulk density, shape, and rolling moment. A total of 4 surrogate munitions and 9 replicas were deployed at each of two water depths adjacent to the quadpod instrument frames (Figure 3). Table 1 shows the complete list of deployed and recovered objects along with brief descriptions and their physical properties.   Divers laid the surrogate and replica munitions on the seafloor around shallow and deep quadpods within the view field of the sector scanning sonar on 21 April 2013 ( Figure 2b). The location and orientation of surrogate and replica munitions were detected by the sector scanning sonar and maintenance diver with video camera. Only objects laid by divers under the shallow quadpod was photographed (Figure 2c). The field of view of the sector scanning sonar is roughly represented by the light blue. The locations of the surrogates are denoted by dark blue circle in the upper left. The other replicas were grouped according to relative bulk density. In this case the red boxes denote the objects that were not recovered from the shallow quadpod site. Thus, the initial surrogate munitions' location and orientation was provided from the TREX13 are only for the shallow quadpod.

Model Description
The open-source Delft3D version 4.04.01 was implemented in the TREX13 area to predict currents, waves, sediment transport, and morphological evolution. Under the wind and tidal forcing, the flow module predicts the currents, feeds the current data into the wave and morphology modules as input, computes the sediment transport, and updates the bathymetry. The flow module [10] uses uniform Chézy-type bottom roughness of 65 m 1/2 s -1 , horizontal eddy viscosity of 0.5 m 2 s -1 , and horizontal eddy diffusivity of 10 m 2 s -1 . The wave module [11] with the stationary mode is used to predict the wave generation, propagation, dissipation, and non-linear wave-wave interactions in the nearshore environment with inputs such as water level, bathymetry, wind, and currents from the flow module. The wave module uses Simulating Waves Nearshore (SWAN), which is a thirdgeneration model derived from the Eulerian wave action balance equation [12]. The coupling time between the flow and wave modules is set to 1 hour. The morphology module works in an integrated way with the wave and flow modules in a cycle. This system is a process-based model that considers the impact of waves, currents, and sediment transport on morphological changes.

Model Grids and Time Steps
Two grids with different grid cell sizes were nested (

Wind and Tidal Forcing
The wind input files were set up using the ERA5 Reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF), with 0.25 (around 28 km) resolution [19] for the flow and wave modules. The Global Inverse Tide Model TPXO 8.0, included in the DDB, was used to create the boundary conditions for the flow module. For the alongshore boundary, the water level with astronomic forcing was imposed. The water level gradient (a so-called Neumann boundary condition) was chosen with a constant zero water level slope in the longshore direction for both across-shore open boundaries. It allows for flow to leave and enter the lateral boundaries with no spurious circulation [20].

Initial and Boundary Conditions
As an initial condition, the water level was set to zero. Additionally, the sediment transport boundary conditions were set by specifying the inflow concentration as zero kg/m 3 . The initial condition for the sand sediment was set as a uniform zero concentration, and the initial bed of sediment was set to 5 m. Wave boundary conditions were set based on the measurements from the deep quadpod location using the significant wave height, wave period, wave directions, and directional spreading. These parameters were applied uniformly on the three open boundaries. The spin-up interval of 720 minutes was established to prevent any influence of a possible initial hydrodynamic instability on the bottom change calculation, which starts only after the spin-up interval. The sediment type was set as sand with a sediment-specific density of 2,650 kg/m 3 .
The calibration was directed to adjust the parameters and allow a better agreement between the model output and measurements. To calibrate the model, water level, waves, and currents from model results were compared with observations. Measurements between 21 and 27 April present a significant variation in water level, waves, and currents. For this reason, this period was selected for calibration. During this process, parameters were adjusted separately. While one was fine-tuned, the others remained constant.

Model Output
The time resolution for the Delft3D output data is 1 hour. The output from the flow module includes the current velocity, Uc = ive+jvn, with (i, j) the unit vectors in longitudinal and latitudinal directions, and the mean water level (h). Let Uc = (ve 2 +vn 2 ) 1/2 be the current speed. The output from the wave module includes the wave peak-period (TP), significant wave height (HS), and wave direction.
The bottom wave orbital velocity (Ubr) results of interactions between surface waves and the seafloor.
A well-established linear wave model with Matlab function [21] is used to calculate the bottom orbital velocity Ubr with the water depth (h), significant wave height (HS), and peak period (Tp) (see Appendix D in [21]). The bottom water velocity vector of combined current and waves is represented by Vw with |Vw| = Uc + Ubr and the orientation, ψ = tan -1 (vn/ve). Figure

Object Mobility Model
Io is the rolling moment of munition about the symmetric axis of the munition (see Figure B1); TF is the forward torque caused by the drag force (Fd) and lift force (Fl) (see Appendix C); pB = B/D, is the percentage burial, and θopb is the cylindrical object burial percentage Shields parameter with large aspect ratio (see Appendix B); (ρo, ρw) are densities of object and water; Π is the volume of the munition. Let the relative horizontal velocity of the rolling object be defined by Substitution of (1) into (2) and use of (3) lead to a special Riccati equation, The special Riccati equation (4) has an analytical solution from integration from tk to tk+1 (k = 0, 1, 2,…, K-1)    (6) with αk and βk as known constants during the integration. Substitution of (6) into (3) Integration of l with respect to time t leads to the munition's displacement.

Object Scour Model
Existing studies on scour burial were all concentrated on motionless objects. The ratio between the fluid force (bottom shear stress) and the weight of the sediment particles, i.e., the sediment Shields parameter (θsed), is crucial for scour burial of motionless object and in turn for prediction of the percentage burial parameter pB(t) = B/D [15]- [16]. Here, f is the wave friction factor [23], ρs the sediment grain density, and d50 the medium sand  As pointed in [23], the equilibrium percentage burial pB,eq for motionless cylinders induced by scour tends to increase as θsed increases. An empirical formula has been established, with different choices of the coefficients (a1, a2, a3) determined experimentally for cylinders subject to steady currents: a1 = 11, a2= 0.5, a3= 1.73 [24], a1= 0.7, a2 = a3 = 0 [25], a1 = 2, a2 = 0.8, a3 = 0 [26], and for cylinders under waves (depending on wave period): a1 = 1.6, a2 = 0.85, a3 = 0 for Tp longer than 4 s [27]. For motionless cylinders before scour burial reaching equilibrium the percentage burial follows an exponential relationship where the e-folding time scale T* is given by With the bottom wave orbital velocity (Ubr), sediment density(ρs), medium grain size (d50), and in turn the sediment Shields parameter (θsed), the equilibrium object percentage burial (pB,eq) are calculated using (9) with coefficients a1 = 1.6, a2 = 0.85, a3 = 0. The sediment supporting depth b (or pb) is calculated from burial depth B (or pB) using (10), i.e., It is noted that the predicted burial percentage (pB) computed from (10) represents the depth that an object on the surface would bury to at that moment. But an object deployed at the beginning of the time sequence would always remain buried at the deepest burial it has reached so far. The burial depth of the base of the object below the ambient seabed is equivalent to the greatest depth that the scour pit has reached up to that point in time [28]. In other words, scour only acts to bury an object deeper. It can never unbury (re-exposure) as the time series suggested by (10). Similar to [28], a simple parameterization was proposed [9] to represent the re-exposure process starting from k (= 1, 2, …): (a) doing nothing if with w the weight coefficient. In this study, we take w = 0.80. The object mobility and burial model consists of Eqs.(6)-(13).

Prediction of Object's Mobility and Burial
Munition's mobility and burial were predicted using the object mobility and burial models with the environmental variables predicted by the Delft3D (Figure 1) as the forcing term. The model was integrated for each surrogate (or replica) deployed in the shallow quadpod with its initial location and orientation (Figure 2c (5)], the object's displacement at the next time step l(tk+1) can be predicted using (7).
Based on the known initial locations of the objects at the shallow quadpod (Figure 2c), the model predicts the object's burial percentage [pB(tk)] shown in Figure 8, the object's percentage burial Shields parameter [θopb(tk)] shown in Figure 9, and the object's displacement [l(tk)] shown in Figure 10. The burial percentages pB for all the objects were less than 0.5 except during the storm event on 12:00 5 May to 00:00 6 May 2013 local time marked by two dotted lines (Figure 8). The red color in Figure 9 shows that the object 's rolling condition [θopb>1] is satisfied. The surrogate and replica munitions' mobility and burial were observed by divers and sector scanning sonar images during the field experiment depicted in Section 2 and in [13].  Figure 11d. Similarity between the observation (Figures. 12a, b) and the model prediction ( Figures. 12 c, d) shows model capability. However, the model-data discrepancy shows the model limitation. For example, yaw of munitions D3 and D6 was observed (Figures 11a, b), but not predicted (Figures 11c, d). The munition D3 (rightmost triangles in  where (eh, ev) are horizontal and vertical unit vectors (see Figure A1). The sediment compressive normal stress FS is decomposed as, With b as the axis of rotation, the sediment above (below) the depth b generates torque to resist (enhance) the rolling with the total torque from the sediment, If we assume that at the depth b the total torque from the sediment is zero (i.e., zero-sum sediment torque for rolling), Here, we take λ = 0.453 in this study.

Appendix B Dynamics of Rolling Object
The drag force (Fd) and lift force (Fl) (see Appendix C) roll the object forward with the torque TF ( Figure B1),   Here, θopb is the object percentage burial Shields parameter based on pb [14]; and So is the relative density of the object. For motionless munition (uo = 0), the condition for the object to move is obtained through substituting (B4) into (B3) where Io is the rolling moment about the symmetric axis of the munition; IA is the rolling moment of munition about the point b (see Figure B1) using the parallel axis theorem; Π is the volume of the munition.

Appendix C Drag, Lift, Buoyancy Forces, and Added Mass
The drag force (Fd), lift force (Fl), buoyancy force (Fw), and added mass (Fa) exerted on the object for rolling by the perpendicular component, U, are given by The vortex shedding from objects is neglected. Besides, the lift coefficient is less certain, we assume where γ is the ratio of lift coefficient versus drag coefficient with γ being taken as value of 0.2.

Appendix D Drag Coefficient
For cylindrical objects, the drag force is decomposed into along and cross axis components. The drag coefficient across-cylinder 's main axis Cd depends on the Reynolds number, where  = 0.8×10 -6 m 2 /s, is the sea water kinematic viscosity; U is the horizontal water velocity perpendicular to the cylinder's main axis; and D is the cylinder's diameter (see Figure 6).
where η=L/D, is the cylinder's aspect ratio.