*Article* **A New 3D Marine Controlled-Source Electromagnetic Modeling Algorithm Based on Two New Types of Quadratic Edge Elements**

**Bin Xiong 1,2, Yuguo Lu 1,2, Hanbo Chen 1,2,\*, Ziyu Cheng 1, Hao Liu <sup>1</sup> and Yu Han <sup>1</sup>**


**Abstract:** This study proposes a new algorithm for a higher-order vector finite element method based on two new types of second-order edge elements to solve the electromagnetic field diffusion problem in a 3D anisotropic medium. To avoid source singularity in the quasistatic variant of Maxwell's function, a secondary field formulation was adopted. The modeling domain was discretized using two types of quadratic edge hexahedral elements, which were obtained using the edge unification method to reduce variables on each side of two conventional quadratic edge elements. Compared with the traditional quadratic element, the number of unknowns that needed to be solved was significantly reduced. The sparse linear equation of the finite element system was solved using an open-source direct solver called MUMPS. The numerical results demonstrated that the proposed algorithm has the same level of accuracy as the conventional vector finite element method and has a significant advantage over it in terms of computational cost.

**Keywords:** vector finite element method; marine controlled-source electromagnetic method; quadratic edge element; hexahedral element; anisotropic

**Citation:** Xiong, B.; Lu, Y.; Chen, H.; Cheng, Z.; Liu, H.; Han, Y. A New 3D Marine Controlled-Source Electromagnetic Modeling Algorithm Based on Two New Types of Quadratic Edge Elements. *Appl. Sci.* **2023**, *13*, 2821. https://doi.org/ 10.3390/app13052821

Academic Editors: Victor Franco Correia and Marek Krawczuk

Received: 9 December 2022 Revised: 29 January 2023 Accepted: 17 February 2023 Published: 22 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The marine controlled-source electromagnetic method (MCSEM) uses controllable artificial field sources to transmit electromagnetic signals in the ocean and measures the electric or magnetic fields at a location far away from the field source to detect the electrical distribution of the media below the seabed. The resistance characteristics of oil and gas reservoirs are the most important properties, which can generate characteristic surface electromagnetic signals. In other words, MCSEM technology can distinguish underground oil and gas from other fluids. In theory, CSEM measurement is to record data with multiple power-receiver offsets, several different frequencies, and a certain power-receiver arrangement. The marine controlled-source electromagnetic (MCSEM) method has become a popular geophysical exploration tool for offshore hydrocarbon (HC) exploration [1–3]. During the last decade, marine CSEM has been widely used to reduce ambiguities in data interpretation and reduce the risk of exploration. With the rapid development of MCSEM exploration instruments, large-scale and high-quality data have been obtained, resulting in the need to develop high-precision MCSEM three-dimensional (3D) data interpretation tools. The 3D inversion algorithm technology of MCSEM is commonly used for data interpretation. Forward modeling is the basis of inversion, and high-precision forward modeling technology is key to achieving high-precision inversion.

The integral equation (IE) [4–6], finite difference (FD) [7–12], and finite element (FE) methods [13–17] are popular numerical techniques for electromagnetic (EM) field modeling. In the FE method, it is easier to use an unstructured mesh to create an irregular body, locally refine the interest region, and coarsen the boundary area of the model domain [15]. With this processing, numerical modeling becomes efficient and effective. Therefore, these

advantages make the FE method more suitable for simulating the complex EM response of irregular models than FD and IE methods.

Compared with the conventional finite element method, which is based on a node, the edge-based finite element method has the advantage that divergence-free conditions are automatically satisfied by appropriately selected basic functions [18]. In recent years, the vector finite element method based on edge elements has been used to solve the electromagnetic field problem, and is known as Maxwell's equations [19,20]. It should be noted that to date, most formulations for the 3D CSEM problem have implemented the lowest order. In addition, some studies have used higher-order elements to solve three-dimensional forward modeling of CSEM, but these were all aimed at tetrahedral elements [14,21]. To the best of our knowledge, no high-order hexahedron element is currently used for threedimensional forward modeling of the controlled-source electromagnetic method. As we know, high-order edge elements can provide accurate solutions but result in an increase in unknown variables and lead to larger computing costs. Fortunately, we can improve the accuracy of the solution at a small computational cost by eliminating some variables when using higher-order elements [22,23].

In this study, we applied the unification method to the edge to reduce one edge variable on each side of the conventional serendipity and the Lagrange-type quadratic hexahedral element, resulting in two new types of second order edge elements. We also introduced a new vector FE method based on the proposed elements to solve the 3D marine CSEM modeling problem in an anisotropic medium. For a complex bathymetry simulation, we used a general hexahedral element to discretize the modeling domain. We validated our code by using several models in isotropic and anisotropic media.

#### **2. Theory**

In geophysical applications, the displacement current is typically ignored when the frequency of the electromagnetic field is low. Maxwell's equations can then be simplified as follows [24]:

$$
\nabla \times \mathbf{E} = \mathbf{i}\omega\mu\_0 \mathbf{H} \tag{1}
$$

$$
\nabla \times \mathbf{H} = \sigma \mathbf{E} + \mathbf{J}\_s \tag{2}
$$

where the space magnetic permeability is denoted by *μ*0, angular frequency is denoted by *ω*, source current is denoted by *J*s, and conductivity tensor is denoted by *σ*:

$$
\sigma = \begin{pmatrix} \sigma\_x & 0 & 0 \\ 0 & \sigma\_y & 0 \\ 0 & 0 & \sigma\_z \end{pmatrix} \tag{3}
$$

where *σx*, *σ<sup>y</sup>* and *σ<sup>z</sup>* are principal conductivities. In this study, we discuss only the principal axis anisotropy case. However, our formulation also works in the case of general anisotropy.

Calculate the curl at the left and right ends of Equation (1), and use Equation (2) to eliminate the magnetic field:

$$
\nabla \times \nabla \times \mathbf{E} - i\omega\mu\_0 \sigma \mathbf{E} = i\omega\mu\_0 \mathbf{J}\_\mathbf{s} \tag{4}
$$

The total and secondary field formulations are two types of formulations that are typically used for solving the 3D MCSEM modeling problem. In contrast to the total field formulation, the secondary field formulation's source term comprises a primary electric field. This type of source term is smoother than the current source used in the total-field formulation. Therefore, in this study, we use a secondary field formulation to address the marine CSEM modeling problem in an anisotropic medium. In this formulation, the total field can be decomposed into two parts: the background (primary field, *E*p) and the anomaly fields (secondary field, *E*s):

$$E = E\_\sf s + E\_\sf p \tag{5}$$

Similarly, the conductivity tensor can be decomposed into background (*σ*b) and anomalous conductivities (Δ*σ*):

$$
\sigma = \sigma\_{\mathfrak{b}} + \Delta \sigma \tag{6}
$$

With substitution of Equations (5) and (6) into Equation (4), the Helmholtz equation with the electric component of the secondary field can be written as:

$$
\nabla \times \nabla \times \mathbf{E}\_{\sf s} - i\omega\mu\_0 \sigma \mathbf{E}\_{\sf s} = i\omega\mu\_0 \Delta \sigma \mathbf{E}\_{\sf P} \tag{7}
$$

The secondary magnetic component can easily be calculated using Equation (8) once the secondary electric component in Equation (7) can be solved:

$$H\_{\rm s} = \frac{1}{i\omega\mu\_0} \nabla \times E\_{\rm s} \tag{8}$$

#### **3. Edge-Based Finite Element Analysis**

In this study, we used the edge-based FE method with quadratic elements, as discussed in the previous section. We used the edge shape function *N<sup>i</sup>* to approximate the secondary electric field within an edge element:

$$E\_s = \sum\_{i}^{N\_{\text{odge}}} N\_i E\_{s,j} \tag{9}$$

where *N*edge denotes the total number of edges in the edge element. For the first-order edge element, the *N*edge = 12; for conventional serendipity and Lagrange-type second-order hexahedral elements, *N*edge = 36 and 54, respectively.

By substituting (9) into (7) and applying the Galerkin method, we can determine that the weak form of the original differential equation is:

$$\mathbf{R}\_{i} = \int\_{\Omega} \mathbf{N}\_{i} \cdot \left[ \nabla \times \nabla \times \mathbf{E}\_{\sf s} - i\omega\mu\_{0}\sigma \mathbf{E}\_{\sf s} - i\omega\mu\_{0}\Delta\sigma \mathbf{E}\_{\sf P} \right] \mathrm{d}v \tag{10}$$

where Ω denotes the modeling area.

Applying first vector Green's theorem to Equation (10), the secondary electric field is given, and it is continuous across the inner boundary; the surface term in Equation (10) vanishes. Then, the discretized form of Equation (10) for each element can be expressed as follows:

$$\mathbf{R}\_{i}^{\rm N\varepsilon} = \sum\_{1}^{\rm N\varepsilon} \left[ \mathbf{K}^{\rm c} E\_{\rm s}^{\rm c} - i\omega\mu\_{0} \mathbf{M}\_{1}^{\rm c} E\_{\rm s}^{\rm c} - i\omega\mu\_{0} \mathbf{M}\_{2}^{\rm c} E\_{\rm F}^{\rm c} \right] \tag{11}$$

where *K<sup>e</sup>* and *M<sup>e</sup>* <sup>1</sup> and *<sup>M</sup><sup>e</sup>* <sup>2</sup> are the local stiffness matrices defined as follows:

$$\mathbf{K}\_{i,j}^{\epsilon} = \int\_{\Omega\_{\mathbf{f}}} (\nabla \times \mathbf{N}\_i^{\mathbf{e}}) \cdot \left(\nabla \times \mathbf{N}\_j^{\mathbf{e}}\right) \mathrm{d}v \tag{12}$$

$$\mathbf{M}\_{1i,j}^{\mathbf{e}} = \int\_{\Omega\_v} \mathbf{N}\_i^{\mathbf{e}} \cdot \boldsymbol{\sigma} \cdot \mathbf{N}\_j^{\mathbf{e}} \mathrm{d}v \tag{13}$$

$$\mathbf{M}\_{2i,j}^{\mathbf{e}} = \int\_{\Omega\_{\mathbf{e}}} \mathbf{N}\_{i}^{\mathbf{e}} \cdot \boldsymbol{\Delta} \boldsymbol{\sigma} \cdot \mathbf{N}\_{j}^{\mathbf{e}} \mathrm{d}v \tag{14}$$

where Ω<sup>e</sup> denotes the domain of an element.

By assembling a global matrix with local matrices, a linear equation for the elements is obtained. The equation is typically sparse:

$$A\mathfrak{e} = \mathfrak{b},\tag{15}$$

Here *A* denotes the global stiffness matrix, and *b* denotes the RHSs of the equations.

To solve Equation (14), a proper condition should be added to the equation. For simplicity, let us consider the homogeneous Dirichlet boundary condition:

$$\mathfrak{e}|\_{\Gamma} = 0 \tag{16}$$

Here, Γ is the boundary area.

IDR(s) iterative solvers with an ILU preconditioner were used to solve the linear equations in Equation (15).

#### **4. First and Conventional Quadratic Edge Elements**

In the edge-based FE method, we can use regular rectangular, general tetrahedral, general hexahedral, or other complex elements to discretize the modeling domain. To simplify the problem and simulate bathymetry, hexahedral elements were used to discretize the modeling domain in this study.

To calculate the stiffness matrix of the hexahedral element, a general hexahedral element should be transformed into an origin-centered cubic element. A general hexahedron in Cartesian coordinates *xyz* is shown in Figure 1a, and a transformed reference domain (cubic element) in Cartesian coordinates *ξζη* is shown in Figure 1b.

The transformation can be described by the following formulas [25]:

$$\mathbf{x} = \sum\_{i=1}^{8} N\_i^{\varepsilon}(\xi, \eta, \zeta) x\_i^{\varepsilon} \tag{17}$$

$$y = \sum\_{i=1}^{8} N\_i^{\varepsilon}(\mathbb{\zeta}, \eta, \mathbb{\zeta}) y\_i^{\varepsilon} \tag{18}$$

$$z = \sum\_{i=1}^{8} N\_i^c(\mathbb{f}, \eta, \mathbb{f}) z\_i^c \tag{19}$$

where *N<sup>e</sup> <sup>i</sup>* is the scalar node-based shape function, which can be defined as follows:

$$N\_i^{\mathbb{f}}(\xi, \eta, \zeta) = \frac{1}{8} (1 + \xi \mathbb{f}\_i)(1 + \eta \eta\_i)(1 + \zeta \mathbb{f}\_i) \tag{20}$$

where *i* is a local node index of the element.

In the analysis of the vector finite element, two volume integrals (12)–(14) must be evaluated. For a hexahedron element, it is easy to transform with two integrals in the *ξηζ*-coordinate using the Jacobian transform:

$$\mathcal{K}\_{i,j}^{\varepsilon} = \int\_{-1}^{1} \int\_{-1}^{1} \int\_{-1}^{1} (\nabla \times \mathcal{N}\_{i}^{\varepsilon}(\underline{\xi}, \eta\_{\prime} \mathcal{J}\_{\prime})) \cdot \Big(\nabla \times \mathcal{N}\_{j}^{\varepsilon}(\underline{\xi}, \eta\_{\prime} \mathcal{J}\_{\prime})\Big) |I(\underline{\xi}, \eta\_{\prime} \mathcal{J}\_{\prime})| \, \mathrm{d}\underline{\xi} \, \mathrm{d}\eta \, \mathrm{d}\underline{\zeta},\tag{21}$$

$$M\_{1i,j}^{\epsilon} = \int\_{-1}^{1} \int\_{-1}^{1} \int\_{-1}^{1} \mathbf{N}\_i^{\epsilon}(\underline{\zeta}, \eta\_{\prime}\underline{\zeta}) \cdot \boldsymbol{\sigma} \cdot \mathbf{N}\_j^{\epsilon}(\underline{\zeta}, \eta\_{\prime}\underline{\zeta}) |\mathbf{J}(\underline{\zeta}, \eta\_{\prime}\underline{\zeta})| \, \mathrm{d}\underline{\zeta} \, \mathrm{d}\eta \, \mathrm{d}\underline{\zeta},\tag{22}$$

$$M\_{2\mathbf{i},\mathbf{j}}^{\varepsilon} = \int\_{-1}^{1} \int\_{-1}^{1} \int\_{-1}^{1} \mathbf{N}\_{i}^{\varepsilon}(\underline{\xi}, \eta, \underline{\zeta}) \cdot \Delta \boldsymbol{\sigma} \cdot \mathbf{N}\_{\underline{j}}^{\varepsilon}(\underline{\xi}, \eta, \underline{\zeta}) |I(\underline{\xi}, \eta, \underline{\zeta})| \, \mathrm{d}\underline{\xi} \, \mathrm{d}\eta \, \mathrm{d}\underline{\zeta},\tag{23}$$

*J*(*ξ*, *η*, *ζ*) and |*J*(*ξ*, *η*, *ζ*)| denote the Jacobian matrix and determinant of*J*(*ξ*, *η*, *ζ*), respectively.

$$J(\xi, \eta, \zeta) = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta} \\ \frac{\partial x}{\partial \zeta} & \frac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{bmatrix} \tag{24}$$

In the next section, we introduce several conventional first-order and quadratic hexahedral elements and two quadratic edge hexahedral elements of the new type.

#### *4.1. First-Order Edge Hexahedral Element*

A conventional first-order hexahedral element is shown in Figure 2, which has eight nodes and twelve edges. The shape function of the edge on the sides along the *ξ*-axis direction are accordingly defined as [25]:

$$N\_i = \frac{1}{8}(1 + \eta\_i \eta)(1 + \zeta\_i \zeta)\Delta\xi \tag{25}$$

**Figure 2.** Conventional first-order edge hexahedral element.

Similarly, the shape function of the edge on the sides along the *η*-axis or *ζ*-axis direction is defined as:

$$N\_{\vec{j}} = \frac{1}{8} \left( 1 + \xi\_{\vec{j}} \xi \right) \left( 1 + \zeta\_{\vec{j}} \zeta \right) \Delta \eta \tag{26}$$

and

$$\mathbf{N}\_k = \frac{1}{8} (1 + \tilde{\varsigma}\_k \tilde{\varsigma}) (1 + \tilde{\varsigma}\_k \tilde{\varsigma}) \Delta \tilde{\varsigma} \tag{27}$$

where *ξ*, *η* and *ζ* denote the local coordinate in the elements.

#### *4.2. A Quadratic Hexahedron Element of a Conventional Serendipity Type*

Figure 3 shows a conventional serendipity quadratic hexahedron element with 12 edges on its surface, 20 nodes in the element, and 24 edges on the sides. The shape function for the serendipitous quadratic hexahedron element can be defined as follows: the shape function of the edge on sides which along *ξ*-axis direction should be:

*<sup>N</sup><sup>i</sup>* <sup>=</sup> <sup>1</sup> 8 (1 + *ηiη*)(1 + *ζiζ*)(4*ξiξ* + *ηiη* + *ζiζ* − 1)Δ*ξ* (28)

**Figure 3.** The location and number of shape functions of a conventional serendipity quadratic hexahedron element (SC). The green arrows and the red arrows represent the shape functions on the edge and on the face, respectively.

The edge shape function of the edges along the *ξ*-axis direction on the *ζ* = ±1 surface:

$$N\_i = \frac{1}{4} \left( 1 - \eta^2 \right) \left( 1 + \zeta\_i \zeta\_i \right) \Delta \xi \tag{29}$$

The edge shape function of the edges along the *ξ*-axis direction on the *η* = ±1 surface:

$$N\_i = \frac{1}{4} \left( 1 - \zeta^2 \right) (1 + \eta\_i \eta) \Delta \xi \tag{30}$$

The sides along the *η*-axis direction on the sides:

$$N\_{\circ} = \frac{1}{8} \left( 1 + \xi\_{\circ}\xi \right) \left( 1 + \xi\_{\circ}\xi \right) \left( 4\eta\_{\circ}\eta + \xi\_{\circ}\xi + \xi\_{\circ}\xi - 1 \right) \Delta\eta \tag{31}$$

The edges along the *η*-axis direction on the *ξ* = ±1 surface:

$$N\_{\circ} = \frac{1}{4} \left( 1 - \zeta^2 \right) \left( 1 + \xi\_{\circ} \xi\_{\circ} \right) \Delta \eta \tag{32}$$

The edges along the *η*-axis direction on the *ζ* = ±1 surface:

$$N\_{\circ} = \frac{1}{4} \left( 1 - \zeta^2 \right) \left( 1 + \zeta\_{\neq} \zeta \right) \Delta \eta \tag{33}$$

The edges along the *ζ*-axis direction on the *η* = ±1 surface:

$$\mathbf{N}\_k = \frac{1}{4} \left( 1 - \zeta^2 \right) (1 + \eta\_k \eta) \Delta \zeta \tag{34}$$

#### *4.3. Lagrange Quadratic Hexahedron Element*

Figure 4 shows a Lagrange quadratic hexahedron element (LC) with 27 nodes and 24 edges on the sides, 24 edges on the surface, and 8 edges within the element. The shape function of the Lagrange quadratic hexahedron element is defined as the edges on the sides along the *ξ*-axis direction:

$$\mathbf{N}\_{i} = \frac{1}{8} (1 + 4\xi\_{i}\xi)\eta\zeta(\eta\_{i} + \eta)(\zeta\_{i} + \zeta)\Delta\xi\tag{35}$$

**Figure 4.** The location and number of shape functions of a Lagrange quadratic hexahedron element (LC). The green arrows and the red arrows represent the shape functions on the edges and on the faces, respectively.

The edges along the *ξ*-axis direction on the *ζ* = ±1 surface:

$$\mathbf{N}\_{i} = \frac{1}{4} (1 + 4\xi\_{i}\xi\_{i})\zeta \left(1 - \eta^{2}\right) (\zeta\_{i} + \zeta) \Delta\xi \tag{36}$$

The edges along the *ξ*-axis direction on the *η* = ±1 surface:

$$\mathbf{N}\_{i} = \frac{1}{4} (1 + 4\xi\_{i}\overline{\mathfrak{z}}) \eta \left(1 - \overline{\mathfrak{z}}^{2}\right) (\eta\_{i} + \eta) \Delta \overline{\mathfrak{z}} \tag{37}$$

The edges along the *ξ*-axis direction within the element:

$$N\_i = \frac{1}{2} (1 + 4\xi\_i\overline{\xi}) \left(1 - \eta^2\right) \left(1 - \overline{\zeta}^2\right) \Delta\overline{\xi} \tag{38}$$

The edges along the *η*-axis direction on the sides:

$$\mathbf{N}\_{\circ} = \frac{1}{8} \left( 1 + 4\eta\_{\circ}\eta \right) \xi \zeta \left( \xi\_{\circ} + \xi \right) \left( \zeta\_{\circ} + \zeta \right) \Delta \eta \tag{39}$$

The edges along the *η*-axis direction on the *ζ* = ±1 surface:

$$N\_{\circ} = \frac{1}{4} (1 + 4\eta\_{\circ}\eta) \zeta \left(1 - \zeta^2\right) \left(\zeta\_{\circ} + \zeta\right) \Delta\eta \tag{40}$$

The edges along the *η*-axis direction on the *ξ* = ±1 surface:

$$N\_{\circ} = \frac{1}{4} \left( 1 + 4\eta\_{\circ}\eta \right) \xi \left( 1 - \zeta^2 \right) \left( \xi\_{\circ} + \xi \right) \Delta\eta \tag{41}$$

The edges along the *η*-axis direction within the element:

$$N\_{\circ} = \frac{1}{2} \left( 1 + 4\eta\_{\circ}\eta \right) \left( 1 - \xi^2 \right) \left( 1 - \xi^2 \right) \Delta\eta \tag{42}$$

The edges along the *ζ*-axis direction on the sides:

$$N\_k = \frac{1}{8} (1 + 4\zeta\_k \zeta) \zeta \eta \left(\zeta\_k + \zeta\right) (\eta\_k + \eta) \Delta \zeta \tag{43}$$

The edges along the *ζ*-axis direction on the *ξ* = ±1 surface:

$$\mathbf{N}\_{k} = \frac{1}{4} (1 + 4\zeta\_{k}\zeta)\zeta \left(1 - \eta^{2}\right) (\xi\_{k} + \xi) \Delta \zeta \tag{44}$$

The edges along the *ζ*-axis direction on the *η* = ±1 surface:

$$\mathbf{N}\_{k} = \frac{1}{4} (1 + 4\zeta\_{k}\zeta)\eta \left(1 - \zeta^{2}\right) (\eta\_{k} + \eta)\Delta\zeta \tag{45}$$

The edges along the *ζ*-axis direction within the element:

$$\mathbf{N}\_k = \frac{1}{2} (1 + 4\zeta\_k \zeta) \left(1 - \eta^2\right) \left(1 - \xi^2\right) \Delta \zeta \tag{46}$$

#### *4.4. A New Type of Quadratic Edge Element*

There are two edges on the side of a conventional quadratic hexahedron element. Within an edge element, because the curls of the edge shape functions are linearly dependent, many edge variables are redundant. The computation time and storage requirements of these elements has increased significantly. Therefore, eliminating these redundant variables helps improve efficiency at no cost to computational accuracy. In this section, we adopt a new method to eliminate redundant variables from the conventional quadratic element.

Let us first consider a one-sided conventional quadratic edge element, as shown in Figures 3 and 4. We named the two edges on this side edge1 and edge2. *N*<sup>1</sup> and *N*<sup>2</sup> are the shape functions of edge1 and edge2, and their corresponding edge variables are *A*<sup>1</sup> and *A*2. Because one of these is redundant, we decided to eliminate *A*<sup>2</sup> and consider *A*<sup>1</sup> = *A*2. Assume *A*side is the vector field interpolation contribution of the two edges:

$$\begin{array}{l} A\_{\text{side}} = \mathbf{N}\_1 A\_1 + \mathbf{N}\_2 A\_2 \\ = (\mathbf{N}\_1 + \mathbf{N}\_2) A\_1 \end{array} \tag{47}$$

Make a line integration of *A* along this side and name it as *A*unif. Note that the following integration is orthogonal:

$$\int\_{\varepsilon KL} \mathbf{N}\_{IL} \, \mathbf{ds} = \delta\_{(IL)(KL)} \tag{48}$$

Here, d*s* is the differential vector of the line element along edge *eKL*, and *δ*(*IL*)(*KL*) represents the Kronecker delta. Thus, the integration *A*unif is written as:

$$A\_{\rm surf} = \int\_{\rm edge1} A \cdot \mathbf{ds} + \int\_{\rm edge2} A \cdot \mathbf{ds} = A\_1 + A\_2 = 2A\_1 \tag{49}$$

Then Equation (47) is rewritten as:

$$A\_{\rm side} = (\mathbf{N}\_1 + \mathbf{N}\_2)A\_1 = (\mathbf{N}\_1 + \mathbf{N}\_2)\frac{A\_{\rm unif}}{2} = \mathbf{N}\_{\rm unif}A\_{\rm unif} \tag{50}$$

In Equation (50), we obtain a new shape function *N*unif = (*N*<sup>1</sup> + *N*2)/2 for the edge, as shown in Figure 5b. According to this, *N*unif can be written in detail as follows:

**Figure 5.** Two-edge unification: (**a**) edge before unification; (**b**) edge after unification.

For the serendipity quadratic hexahedron element, the edges on the sides along the *ξ*-axis direction:

$$N\_k = \frac{1}{2} (1 + 4\zeta\_k \zeta) \left(1 - \eta^2\right) \left(1 - \zeta^2\right) \Delta \zeta \tag{51}$$

The edges on the sides along the *η*-axis direction:

$$N\_{\circ} = \frac{1}{8} \left( 1 + \zeta\_{\circ}\zeta\right) \left( 1 + \zeta\_{\circ}\zeta\right) \left( \zeta\_{\circ}\zeta + \zeta\_{\circ}\zeta - 1 \right) \Delta\eta \tag{52}$$

The edges on the sides along the *ζ*-axis direction:

$$N\_k = \frac{1}{8}(1 + \xi\_k \xi)(1 + \eta\_k \eta)(\xi\_k \xi + \eta\_k \eta - 1)\Delta\xi \tag{53}$$

For the Lagrange quadratic hexahedron element, the edges on the sides along the *ξ*-axis direction:

$$\mathcal{N}\_i = \frac{1}{8} \eta \mathcal{J} (\eta\_i + \eta)(\mathcal{J}\_i + \mathcal{J}) \Delta \xi \tag{54}$$

The edges on the sides along the *η*-axis direction:

$$\mathbf{N}\_{\circ} = \frac{1}{8} \xi \zeta \left(\xi\_{\circ} + \xi\right) \left(\zeta\_{\circ} + \zeta\right) \Delta \eta \tag{55}$$

The edges on the sides along the *ζ*-axis direction:

$$\mathbf{N}\_{k} = \frac{1}{8}\xi\eta\left(\xi\_{k} + \xi\right)(\eta\_{k} + \eta)\Delta\zeta\tag{56}$$

The shape functions of the edges in the element and the faces remained intact. Therefore, we obtained two new quadratic edge elements (Figure 6). There is only one edge on the sides of the two new quadratic edge elements (Figure 6); therefore, the redundant variable can be effectively eliminated. Table 1 lists various types of elements and the distribution of their basis functions.

**Figure 6.** The location and number of shape functions of two new types quadratic elements: (**a**) new serendipity quadratic element (SN); (**b**) new Lagrange quadratic element (LN). The green arrows and the red arrows represent the shape functions on the edges and on the faces, respectively.

**Table 1.** Elements and the distribution of their basic functions.


#### *4.5. Interpolation and Convergence*

Interpolation capacity and convergence tests are performed using two new types of second-order element. Those functions are approximated in the domain defined by *Ω* = [−1 1] 3 . Figure 7 shows the error of interpolation of the tested functions and their respective convergence rates. The two basis functions (SC, LC) have similar convergence rates, both close to one.

**Figure 7.** Interpolation errors and convergence rates for the two basis functions approximated by EFEM. H is the length of edge of a element.

#### **5. Model Studies**

All MCSEM response simulations were performed on a platform with hardware and operation, as follows. The CPU was a double Intel processor, each processor having 6 cores and 12 threads, and RAM memory was 64 GB at a frequency of approximately 1333 MHz. The operating system was a 64-bit Windows 7 professional version.

#### *5.1. Verification of the New Algorithm*

To verify the new algorithm, an isotropic horizontal layered geoelectrical model, as shown in Figure 8, was considered first. In this model, air and seawater were 1000 m thick. A 100 m thin-layer reservoir was embedded in the sediment at 1500 m to 1600 m depth. The EM field response was based on the model presented in Table 2. An *x*-oriented 1 Hz frequency current source was placed at the location (0, 0, 950 m). The receiving stations were located on the seafloor (*z* = 1000 m). The modeling domain was 3 km × 3 km square and 4 km in the vertical direction (−1 km ≤ z ≤ 3 km). A non-uniform rectangular grid was used to discretize Ω. We refined the grid near the current source, receiving stations, and target layer. The analytical response was calculated using the numerical integration of the Hankel transform [25–29].

**Table 2.** The conductivities of an isotropic horizontally layered geoelectric model.


**Figure 8.** Horizontally layered geoelectric model with rectangular mesh; the red star denotes the electric dipole source. The blue, green, red and gray areas represent the air layer, sea water, oil and gas, and bedrock, respectively.

Figures 9 and 10 show the difference in the *x*-component of the electric field by applying conventional serendipity (SC)-and Lagrange (LC)-type second-order hexahedral elements, and the proposed elements (new serendipity (SN) and Lagrange element (LN)). These new quadratic edge elements were at the same level of accuracy as conventional quadratic elements, both higher than that of linear edge elements.

**Figure 9.** A comparison of the *x*-component of the electric field using different FEM methods (serendipity-type second-order element) with the analytical solution.

**Figure 10.** A comparison of the *x*-component of the secondary electric field by using different FEM methods (Lagrange-type second-order element) with the analytical solution.

Table 3 lists the CPU times for different elements. Less time was required when using new quadratic elements than when using conventional quadratic elements. Therefore, one can say that the computation is reduced at no cost to accuracy.

**Table 3.** The comparison of CPU time by using different methods.


*5.2. Off-Shore Hydrocarbon ReservoirMmodel*

Consider a geoelectric model with a hydrocarbon reservoir embedded in sediment, as shown in Figures 11 and 12. In this model, air and seawater were 1000 m thick, similar to the previous model. The conductivity values of each component are listed in Table 4. The modeling domain was a 4 km × 4 km × 4 km cube with 1 km ≤ *z* ≤ 3 km vertically . The hydrocarbon reservoir was buried 2.0–2.1 km deep and occupied a 1 km × 1 km area horizontally. An *x*-oriented 1 Hz frequency current source was placed at the location (−3000, 0, 950) m. Receiving stations were placed on the seafloor (*z* = 1000 m). We used the same type of rectangular grid used in the previous model to discretize the model domain. We also refined the local grid for the current source, receiving stations, and the target layer area.

**Figure 11.** Horizontally layered geoelectrical model with rectangular mesh. The red star identifies the location of the current source. The dark blue, light blue, yellow, and brown areas represent the air layer, sea water, reservoir, and sediment, respectively.

**Figure 12.** Plane view of the grid. The yellow area is the area occupied with reservoir, and the red star identifies the location of the current source.


**Table 4.** The conductivities of a hydrocarbon reservoir geoelectrical model.

Figures 13 and 14 show the *x*-component of the secondary electric field calculated with different edge elements. The results of the four types of second-order elements were in agreement with each other. Furthermore, a strong anomalous field distortion caused by the background conductivity anisotropy can also be observed. Table 5 shows the cost-time of calculation with four types of quadratic edge elements. Similar to the previous model, the efficiencies of the two types of quadratic edge elements (SN and LN) were higher than those of conventional serendipity or Lagrange-type quadratic hexahedron elements (SC, LC).

**Figure 13.** Secondary electric component in the *x*-direction, calculated by using different elements for isotropic sediment.

**Figure 14.** Secondary electric component in the *x*-direction, calculated by using different elements for anisotropic sediment.


**Table 5.** The cost-time of CPU when using different methods.

#### *5.3. Pyramid Type Bathymetry Model*

Further verification of the new algorithm was performed by considering a pyramidtype bathymetry model (Figure 15) without a reservoir. The top area of the pyramid-type bathymetry model was 0.5 km × 0.5 km in square, at 990 m depth. The bottom area was 1.0 km × 1.0 km square, at 1000 m in depth. Air and sea water were 1000 m thick. In the domain of pyramid-type bathymetry, the depth of seawater ranges from 900 m to 1000 m. Figure 16 shows the pyramid-type bathymetry model with a reservoir. The domain of the reservoir was Ω = {−1000 m ≤ *x*, *y* ≤ 1000 m, 1500 m ≤ *z* ≤ 1600 m}. The modeling domain was Ω = {−4 km ≤ *x*, *y* ≤ 4 km, −1 km ≤ *z* ≤ 3 km}. In these models, the conductivities of air, seawater, and sediment were the same as previously, whereas the reservoir's conductivity was 0.05 s/m. The electric dipole source was *x*-oriented at (−3000, 0, 800) m at 1 Hz frequency, and the receiving station was located at *z* = 800 m.

**Figure 15.** Section view of the grid used for the bathymetry model without reservoir. The red star locates the dipole source. The dark blue, pink, blue, and green areas represent the air layer, sea water, oil and gas, and bedrock, respectively.

**Figure 16.** Section view of the grid used for the bathymetry model with reservoir. The reservoir occupies the blue area, and the current source is located at the red star. The dark blue, pink, blue, and green areas represent the air layer, sea water, oil and gas, and bedrock, respectively.

Figures 17–20 show the secondary electric component in the *x*-direction of the bathymetry model, which has no reservoir or one reservoir in the isotropic sediment. Figures 21–24 show a comparison of the *x*-component of the secondary electric field of the bathymetry model with and without a reservoir in anisotropic sediment. Bathymetry and anisotropic background clearly affected the distribution of the marine controlled-source electromagnetic field (Figures 18–24), which should be considered in the processing of MCSEM data. Figures 25 and 26 show a comparison of the degrees of freedom (DOF) and CPU time using different elements, respectively. Clearly, these two new elements contributed significantly to reducing the CPU time in computing.

**Figure 17.** No-reservoir bathymetry model in isotropic sediment. The secondary electric field was calculated by using different second-order edge elements; here is the comparison.

**Figure 18.** No-reservoir bathymetry model in isotropic sediment. Amplitude comparison between background field (blue line) and total field (dashed red line) for *Ex* component, at *y* = 0 m.

**Figure 19.** Reservoir-embedded bathymetry model in isotropic sediment. Second electric component in the *x*-direction, which was calculated by using different quadratic edge elements.

**Figure 20.** Reservoir-embedded bathymetry model in isotropic sediment. Amplitude comparison between background field (blue line) and total field (dashed red line) for *Ex* component, at *y* = 0 m.

**Figure 21.** No-reservoir bathymetry model in anisotropic sediment. Second electric component in the *x*-direction, which was calculated by using different quadratic edge elements.

**Figure 22.** No-reservoir bathymetry model in anisotropic medium. Amplitude comparison between background field (blue line) and total field (dashed red line) for *Ex* component, at *y* = 0 m.

**Figure 26.** Comparison of CPU time using different elements. Model 1 = bathymetry + isotropic sediment; model 2 = bathymetry + isotropic sediment + reservoir; model 3 = bathymetry + anisotropic sediment; model 4 = bathymetry + anisotropic sediment + reservoir.

#### **6. Conclusions**

We developed a new vector finite element algorithm to solve the 3D marine CSEM modeling problem in an anisotropic medium. An algorithm based on a secondary field formulation and two new quadratic edge elements was proposed. Compared with the vector FEM algorithm based on the conventional quadratic edge element, the proposed algorithm significantly reduced the computational cost without losing accuracy. The algorithm in this paper has good generality and is suitable for solving other electromagnetic induction problems in isotropic and anisotropic media, including airborne and borehole electromagnetic methods.

**Author Contributions:** Conceptualization, B.X. and H.C.; methodology, H.C.; software, H.C. and Y.L.; validation, Z.C. and Y.H.; formal analysis, Y.L.; investigation, H.L.; resources, H.C.; data curation, H.C.; writing—original draft preparation, H.C.; writing—review and editing, Y.L.; visualization, B.X.; supervision, B.X. and H.C.; project administration, B.X.; funding acquisition, H.C. and B.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China No. 42174080. It was also supported by the Natural Science Foundation of Guangxi (No. 2020GXNS-FAA297079) and the China Postdoctoral Science Foundation funded project (No. 2021MD703820) and the Research Start-up Foundation of Guilin University of Technology (No. RD2100002165).

**Institutional Review Board Statement:** The study did not require ethical approval.

**Informed Consent Statement:** The study did not involve humans.

**Data Availability Statement:** The data is unavailable due to privacy or ethical restrictions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Yuxi Wu 1, Kui Zhu 2,\*, Hao Qin 1, Yang Wang 3, Zhaolong Sun 2, Runxiang Jiang 2, Wanhu Wang 3, Jiaji Yi 3, Hongbing Wang <sup>3</sup> and Enjin Zhao <sup>1</sup>**


**Abstract:** The coastline of Shantou is tortuous, while the hydrodynamic environment is complicated. In this paper, the hydrodynamic model is established by the FVCOM (Finite Volume Coastal Ocean Model); the open boundary conditions such as water level, river, and wind field are the input; and the model is verified by tidal harmonic function. According to the previous research, the typhoon wind field with a 10-year return period is selected for storm surge simulation. When there is a bank, the accumulated water on the land cannot enter the ocean due to the block of the bank but accumulates on the inner side of the bank, resulting in higher accumulated water, but less than 0.5 m. In the aspect of sediment deposition, a sediment transport model is established to analyze the sediment deposition in Shantou Port and its surrounding waters. Some reasonable suggestions are put forward for the sediment deposition in Shantou. According to the simulation results, the following conclusions can be drawn: (1) In the case of typhoon storm surge in the return period of 10 years, the bank can effectively protect the inland. Still, accumulated water will collect near the bank. (2) The offshore water level will rise by 0.4 m after adding a bank. (3) The sediment in Shantou Bay mainly comes from the ocean sediment caused by tides, and the largest sedimentation occurs in the main channel.

**Keywords:** numerical simulation; hydrodynamic force; sediment transport; embankment

#### **1. Introduction**

The coastline of China is twisty and turning, and foreign trade in coastal areas prospers. The economy develops rapidly but is also threatened by marine disasters [1], among which the storm surge brings the most severe losses [2–4]. The post-disaster reconstruction after the storm surge requires a large number of materials and personnel investment. Storm surges are generally caused by typhoons [5]. In summer, the southeast coast of China is easily attacked by tropical cyclones, causing massive casualties. As a typical coastal special economic zone in China, Shantou has experienced many typhoons landing [6]. The accompanying storm surges cause heavy financial losses and often cause thousands or even tens of thousands of casualties.

The formation and development mechanism of storm surge was put forward in the early 20th century. Firstly, the theoretical storm surge model was put forward, which was constantly supplemented and improved, and its mechanism was deeply discussed [7]. With the rapid development of geodetic techniques, the research on numerical models has been deepened. Since the 21st century, many numerical models for storm surge models have existed. For example, the POM (Princeton Ocean Model) and the improved ECOM (Estuary Coast and Ocean Model) of Princeton University in the United States are three-dimensional baroclinic models; the Delft-3D model of Delft Water Conservancy Research Institute in the Netherlands; HAMSOM (Hamburg Shelf Ocean Model) developed by the University of Hamburg, Germany; and the FVCOM (An Unstructured Grid Finite Volume Coastal Ocean

**Citation:** Wu, Y.; Zhu, K.; Qin, H.; Wang, Y.; Sun, Z.; Jiang, R.; Wang, W.; Yi, J.; Wang, H.; Zhao, E. Numerical Investigation on the Influence of Breakwater and the Sediment Transport in Shantou Offshore Area. *Appl. Sci.* **2023**, *13*, 3011. https:// doi.org/10.3390/app13053011

Academic Editor: Jianhong Ye

Received: 4 January 2023 Revised: 6 February 2023 Accepted: 15 February 2023 Published: 26 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Model) of the Massachusetts Institute of Technology are all models used previously [8–10]. The storm surge model in this paper is based on the FVCOM model and is coupled with the typhoon data of a 10-year return period has been verified in the reference.

In order to resist the storm surge, many scholars put forward coastal port facilities protection measures and disaster reduction measures [11–13]. Currently, there are mainly two kinds of protection against storm surges: one is to set up a bank, and the other is to build an ecological vegetation coast [14]. The former is a traditional protective measure that can weaken the waves caused by the storm surge. It has a positive effect in weakening the storm surge, but its influence range is limited, and it has little impact on the sea area far away from the shoreline [15–17]; the latter is that blueprint of green seawall recently put forward, and the previous study have shown that salt marsh vegetation can attenuate tidal current velocity, and a model that can be used to study the interaction mechanism between wetland vegetation and hydrodynamic forces during storm surge transit has been established [18]. In this paper, we consider the effectiveness of setting up the embankment to weaken the storm surge under the typhoon storm surge model with a 10-year return period.

No matter whether the embankment is set up or the vegetation is covered, it belongs to human activities, which will destroy the original coastline in this area, leading to the reduction of the sea area and the tidal volume in the bay, thus affecting the sediment carrying capacity of the sea area, accelerating the sediment deposition rate in the bay [19], and changing the coastline shape again, forming the positive feedback of sediment deposition. Therefore, coastal cities are affected by human activities in bays, estuaries, and ports, and sediment deposition often occurs. It is necessary to analyze and evaluate the sediment in the port area to protect the topography and ecological environment of the basin. Tides usually control the deposition, so the sediment model needs to be coupled with the hydrodynamic model for overall analysis [20–23]. Currently, the widely used sediment transport models include the MIKE model, MOHID model, ECOM, FVCOM, Delft3D model, etc. [24–26]. The sediment models at the port are complicated and are affected by three factors: topography, hydrodynamics, and sediment conditions. Aiming at the problem of sediment deposition in Shantou Port, many scholars have conducted detailed and precise research on sediment deposition in Shantou Port [27,28]. Still, the research time is relatively old, and due to the interference of human activities, the water area of Shantou Port is constantly shrinking, and the hydrodynamic conditions are continually changing. Therefore, it is necessary to build a new model of sediment deposition in Shantou Port in time.

In this paper, the hydrodynamic model of Shantou offshore is established, and the model is verified by the tidal harmonic constant. Based on the hydrodynamic model, according to the study of the Shantou storm surge model by Tang Yue Zhao, a 10-year return period typhoon storm surge model is established. The situation of storm surge increasing water and causing disasters in the 10-year return period is analyzed in the circumstance of deposition. At the same time, the sediment transport model is established and verified. The sediment deposition in Shantou Bay and surrounding waters are analyzed using the model, and the causes of sediment deposition are analyzed.

#### **2. Model Introduction and Verification**

#### *2.1. Hydrodynamic Model*

This paper uses FVCOM to simulate the ocean hydrodynamic forces, which is a prediction model for the ocean hydrodynamic environment in offshore areas and estuaries [29]. FVCOM adopts the finite volume method, uses unstructured triangular mesh to solve the original equation discretely, and encrypts the research area locally. On the one hand, the finite volume method can discretely solve the variables of the original equation. On the other hand, it can ensure efficient calculation when the equation is integrated, which fully uses the respective advantages of the finite element method and the finite difference method [30]. In the Cartesian coordinate system, the governing equations of FVCOM

include the continuity equation, momentum conservation equation in three directions, temperature and salinity conservation equation, and density equation.

$$\frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \tag{1}$$

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial y} - fv = -\frac{1}{\rho\_0} \frac{\partial P}{\partial x} + \frac{\partial}{\partial z} (K\_m \frac{\partial u}{\partial z}) + F\_u \tag{2}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} + fu = -\frac{1}{\rho\_0} \frac{\partial P}{\partial y} + \frac{\partial}{\partial z} (K\_m \frac{\partial u}{\partial z}) + F\_v \tag{3}$$

$$\frac{\partial P}{\partial z} = -\rho g \tag{4}$$

$$\frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} = \frac{\partial}{\partial z} (K\_h \frac{\partial T}{\partial z}) + F\_T \tag{5}$$

$$\frac{\partial S}{\partial t} + u \frac{\partial S}{\partial x} + v \frac{\partial S}{\partial y} + w \frac{\partial S}{\partial z} = \frac{\partial}{\partial z} (K\_h \frac{\partial S}{\partial z}) + F\_S \tag{6}$$

$$
\rho = \rho(T, S) \tag{7}
$$

Equation (1) is a continuity equation, where *x*, *y*, and *z* represent the positions in the east-west direction, north-south direction, and vertical direction in the Cartesian coordinate system. *u*, *v*, and *w* represent the velocity components in three directions in the Cartesian coordinate system. Equations (2)–(4) are momentum conservation in the east-west direction, north-south direction, and vertical direction, respectively, and *F* represents the parameter of geostrophic deflection force; *ρ* is the density of local seawater and *ρ*<sup>0</sup> is the density of reference position; *P* is pressure; *Km* is the vertical eddy viscosity coefficient; *G* is the acceleration of gravity; *Fu*, *Fv*, *FT*, and *FS* are diffusion terms of momentum, temperature, and salinity in two horizontal directions, respectively. Equations (5)–(7) are the temperature equation, salinity equation, and density equation, respectively. *T* is temperature, *S* is salinity, and *Kh* is the vertical diffusion coefficient of turbulent kinetic energy.

In the actual simulation process, the natural seabed terrain is often uneven. In order to avoid the false flow and diffusion caused by the rough bottom of the model, it is necessary to transform the σ coordinate vertically to smooth the bottom.

$$
\sigma = \frac{z - \zeta}{H + \overline{\zeta}} = \frac{z - \overline{\zeta}}{D} \tag{8}
$$

Equation (8) is the conversion between the σ coordinate and Cartesian coordinate, where *z* is the specific position of the Cartesian coordinate system or the height of the water surface relative to the still water, *H* is the distance from the reference surface to the seabed, and *D* is the height of the total water column. Under the condition of this coordinate system, the variation range of *σ* is from 0, indicated by seawater, to −1 at the bottom of the ocean, which avoids the unevenness of the bottom in the original Cartesian coordinate system.

#### *2.2. Typhoon Model*

According to the previous research, this paper uses the Holland typhoon model to construct the typhoon wind field with a frequency of 10 years. The model typhoon center pressure is 947 hPa, and the maximum wind speed is 47.26 m/s. The governing equation of the model is as follows:

$$P\_{\mathcal{I}} = P\_{\mathcal{C}} + \left(P\_n - P\_{\mathcal{C}}\right) \exp\left(-\left(R\_{\text{max}}/r\right)^B\right) \tag{9}$$

$$V\_r = \left[\frac{B}{\rho\_a} \left(\frac{R\_{\text{max}}}{r}\right)^B (P\_n - P\_c) \exp\left[-\left(\frac{R\_{\text{max}}}{r}\right)^B\right] + \left(\frac{rf}{2}\right)^2\right]^{0.5} - \frac{rf}{2} \tag{10}$$

Equation (9) is the pressure distribution formula of the Holland typhoon model, where *Pc* is the lowest pressure in the typhoon center and *Pn* is the background pressure, at 1013 hPa; *R*max is the maximum wind speed radius of the typhoon; *r* is the distance from somewhere to the typhoon center; *B* is the scale parameter of the typhoon model. Equation (10) is the wind speed distribution formula of the typhoon model, where the air density (*ρa*) is 1.2 kg/m3; *f* is the parameter of Coriolis force.

There are many formulas to determine the maximum wind speed radius of the typhoon, such as Graham, Tang Jian, Jiang Zhihui, etc. [31–33], and all put forward corresponding formulas. In this study, the procedure proposed by Willoughby [34] is adopted to determine the maximum wind speed radius of a typhoon:

$$R\_{\text{max}} = 51.6 \exp\left(-0.0223 V\_{\text{max}} + 0.0281 \varphi\right) \tag{11}$$

where *ϕ* is the latitude of the place where the typhoon center is located.

The maximum wind speed (*V*max) of the typhoon is constructed according to the relationship between wind speed and air pressure proposed by Atkinson-Hollidy [35], and the empirical formula is given:

$$V\_{\text{max}} = 3.7237 \times (P\_n - P\_c)^{0.6065} \tag{12}$$

#### *2.3. Sediment Model*

This paper uses MIKE software's ST (sediment transport) module to calculate sediment transport. The open boundary of the sediment model is driven by the calculation result of the FVCOM hydrodynamic model, and the wind field data are reanalyzed by ECMWF. MIKE FVCOM and Mike FVCOM use unstructured grids and finite volume methods to solve discrete integrals so they can share the same grid system. The governing equation of the sediment model is as follows:

$$\frac{\partial \overline{\varepsilon}}{\partial t} + \mu \frac{\partial \overline{\varepsilon}}{\partial x} + v \frac{\partial \overline{\varepsilon}}{\partial y} = \frac{1}{h} \frac{\partial}{\partial x} (hD\_x \frac{\partial \overline{\varepsilon}}{\partial x}) + \frac{1}{h} \frac{\partial}{\partial y} (hD\_y \frac{\partial \overline{\varepsilon}}{\partial y}) + \frac{F\_S}{h} \tag{13}$$

Equation (13) is the sediment transport formula of the model, where the average sediment concentration is; *u* and *v* are the average vertical velocities in the x and y directions, respectively; *h* is the water depth there; *Dx* and *Dy* are turbulent diffusion coefficients of sediment in two directions, respectively; *FS* is a function term of sediment source.

Sediment source function *FS* is defined by the following formula:

$$F\_{\sf s} = \begin{cases} E\left(\tau\_b / \tau\_{c\sf t} - 1\right)^n & \tau\_b > \tau\_{c\sf t} \\ 0 & \tau\_{c\sf r} \ge \tau\_b \ge \tau\_{c\sf t} \\ \omega\_s c\_b p\_d & \tau\_b < \tau\_{c\sf t} \end{cases} \tag{14}$$

where *E* the sediment scouring coefficient of the bottom of the river is bed; *ω<sup>s</sup>* is the settling velocity of sediment, which is related to the shape, size, and content of sediment particles. *pd* is the possible settlement probability; *τ<sup>b</sup>* is the shear stress at the bottom of the bed; *τce* is the critical shear stress initiated by surface sediment; *τcd* is the acute shear stress at which suspended sediment begins to settle.

The range of the sediment transport model (Figure 1) selects Shantou Port, including the Rongjiang River and Hanjiang River and its surrounding waters. Based on the hydrodynamic model, rivers are added. The study area is narrowed. Thus, the running time of the model is shortened, and the accuracy of the key research area is improved.

**Figure 1.** Grid of sediment transport model.

#### **3. Verification of the Model**

#### *3.1. Verification of Tidal Harmonic Constant*

The tidal harmonic constant of the interpolated hydrodynamic model is verified to prove the validity of the model. In this paper, three long-term tidal stations in the study area are selected, namely Yunaowan Tidal Station (117◦06 E, 23◦24 N), Shantou Tidal Station (116◦44 E, 23◦20 N), and Haimen Tidal Station (116◦37 E, 23◦11 N). In the model of the previous section, the tidal level data of three tide points are derived. The measured tidal level data and the derived tidal level data are calculated by the T-tide program, and the tidal harmonic constants of the corresponding positions are obtained. The simulation results are compared with the measured results' tidal harmonic constants of the four main tidal components K1, O1, M2, and S2. See Table 1 for the verification results of the model. It can be seen from the table that the amplitude error of other tidal components is less than 4 cm, and the late angle error is less than 12 degrees, except for M2, which is slightly larger. The mistake of the Shantou tide gauge station is somewhat more prominent, influenced by its geographical location. Generally speaking, the tidal harmonic constant error between the simulated and the measured results is within an acceptable range, which means that the model has sure accuracy and credibility.

#### *3.2. Verification of Deposition Rate*

The sediment transport model is verified by calculating the change of sediment thickness in Shantou Port in the simulation results (Figure 2). The simulation is from 1 January 2015 to 31 December 2019. Shantou Port's natural sediment deposition and surrounding waters are simulated without human intervention. Through calculation, the average deposition thickness of Shantou Port in five years is 0.416 m, and the average annual deposition thickness is 8.3 cm. The siltation speed of Shantou Port has been impacted by human activity throughout history, as evidenced by available data [27]. It can be roughly estimated that, if the initial conditions, such as water area and sediment content of Shantou Port, are constant, the average deposition thickness from 2015 to 2019 should be about 9 cm, slightly higher than the model simulation results. We are assuming that the initial conditions, such as water area and sediment content of Shantou Port, are constant at the beginning of the People's Republic of China. The average deposition thickness from 2015 to 2019 should be about 9 cm, slightly higher than the model simulation results. However, considering the

breakwater construction, part of the sediment entering the port will settle at the estuary, resulting in the real sedimentation rate being lower than the estimated value. Therefore, the research model is more accurate.


**Table 1.** Verification of tidal harmonic constant.

**Figure 2.** Changes of the bed bottom in the Shantou Port in five years.

#### **4. Model Area and Data**

The study area is from 116◦22 E west to 117◦48 E in the east and from 22◦68 N in the south to 23◦73 N in the north, including Shantou's offshore and surrounding waters. We used SMS (surface-water modeling system) software to lay unstructured triangular grids. The density of open boundary nodes decreases gradually from land to sea, which can improve computational efficiency while encrypting the critical research areas. The grid resolution in the coastal part of the open boundary is high, and the grid resolution in the main study area near the coastal area of Shantou City is up to 100 m. The research area grid contains 55,790 computing nodes and 108,354 irregular triangular grids (Figure 3). Considering that the research object belongs to the ocean surface, it is divided into five *σ* layers vertically.

**Figure 3.** Grid setting in Shantou offshore (the red box in the figure is the coastal area of Shantou City where the grid resolution is up to 100 m).

High-precision data are an essential condition for the model to obtain reliable results. The water depth data in this paper adopt the global seabed topographic data of ETOPO1 provided by the NOAA (National Ocean and Atmospheric Administration) and the chartsounding data with the number 81001. The ETOPO1 data have the highest accuracy among ETOPO data, with a spatial resolution of 1/60◦ × 1/60◦. The bathymetric data of the chart are the measured data, which have a high density in the coastal areas. The datum plane of ETOPO1 data is the global mean sea level, while the datum plane of the chart is the local mean tide height datum plane. After conversion, the tide level starting surface of Shantou City is uniformly used as the datum plane, Figure 4 shows the water level in Shantou offshore and Figure 5 shows the study area after inserting topographic data.

**Figure 4.** Bathymetric chart (The figure shows the water level in Shantou offshore. The colder color means a deeper water level).

**Figure 5.** Terrain of the study area.

#### **5. Results of Typhoon Storm Surge Simulation**

*5.1. Analysis of Typhoon Storm Surge Disasters with Different Intensities*

According to the statistics of tropical cyclones in Shantou since the founding of the People's Republic of China, Tang Yue Zhao designed storm surge models in Shantou with different return periods combined with the Holland typhoon model and simulated storm surge with parametric typhoon models. It was found that southeast-west typhoons were the most common, and the typhoon center pressures corresponding to the return periods of 10 years, 50 years, 100 years, and 1000 years were 947 hPa, 937 hPa, 920 hPa, and 905 hPa, respectively [36]. The maximum water increase in different return periods is obtained from the four stations of Nan'ao Station, Chaoyang District, Chenghai District, and Shantou Port as shown in Figure 6.

**Figure 6.** The maximum water increase in different return periods is obtained from the four sta-tions of Nan'ao Station, Chaoyang District, Chenghai District, and Shantou Bay (The maximum water increase is shows an increasing trend with the increase of return period years).

In actual working conditions, the probability of a storm surge of more than 100 years is low, and the port project cannot cope with the storm surge of a 100-year return period [37,38]. According to the research results of Tang Yue, this paper selects the typhoon moving from southeast to northwest in the 10-year return period and analyzes the protection effect of setting the bank dike on Shantou. To make the simulation results of typhoons and storms more intuitive, the north latitude range is 23◦33 N to 23◦45 N and the west longitude from 116◦66 E to 116◦84 E is selected for local simulation. Some grids in the main urban

area of Shantou extend inland based on the original coastline. At the same time, Shantou's coast is artificial, so the actual shoreline is changed to a bank embankment with a height of 1.5 m. When the external water level is less than 1.5 m, it is regarded as a solid boundary where no water exchange occurs. The standard simulation calculation is carried out when it is higher than 1.5 m. The extended part is shown in Figure 7. The grid extending to the land is in the red line in the figure. The land elevation data are obtained by combining the extracted remote sensing data with the measured elevation data. The red line in the figure is the original shoreline, the dam is set in the model, and its height is set to be 1.5 m higher than the average sea level.

**Figure 7.** Urban coastal grid expansion (The red line in the figure is the original shoreline, the dam is set in the model, and its height is set to be 1.5 m higher than the average sea level).

Figure 8 shows the distribution of water level and flow field in the Shantou sea area and the inundation of storm surge on the coast and land under the southeast-northwest route during the 10-year return period of typhoon storm surge. The center of the typhoon landed in Chenghai District, Shantou City. Under this path, the hurricane hit Shantou City head-on, which has an excellent representative significance. The simulation time is 72 h, of which 48 to 72 h is the transit time of typhoons.

(**b**) The inundation range and depth of Rongjiang River Estuary and main urban area at 2 h during storm surge

(**c**) Water level distribution with velocity vector in the whole study area at 8 h during storm surge

 (**e**) Water level distribution with velocity vector in the whole study area at 14 h during storm surge

(**d**) The inundation range and depth of Rongjiang River Estuary and main urban area at 8 h during storm surge

(**f**) The inundation range and depth of Rongjiang River Estuary and main urban area at 14 h during storm surge

(**h**) The inundation range and depth of Rongjiang River Estuary and main urban area at 20 h during storm surge

**Figure 8.** Tidal current field and inundation of typhoon storm in 10-year return period (The red box on the left is the simulation area on the right).

In Figure 8, the left part is the water level distribution (based on the mean sea level of the Yellow Sea) and velocity vector diagram in the study area, and the right part is the inundation range and inundation depth near the Rongjiang Estuary and the main urban areas. The typhoon reached its maximum inundation 14 h ago, and the deepest accumulated water reached 0.38 m. During the 2 h of storm surge, the typhoon center did not match the study area. Still, under the wind, the water level in the study area obviously increased, especially in the northeast of Nan'ao Island, where the water level rose by more than 0.5 m. At 8 h, the typhoon center entered the study area. The water level gradient in the study area was high near the shore and low at the far shore, and the water level near the Rongjiang Estuary increased by nearly 1 m. At 14 h, the typhoon landed, and the hurricane was located northwest of the study area. The water level was high near the shore, low on the far shore, high in the northeast, and low in the southwest. The submerged area inside the dam expanded, but the accumulated water depth was still within 0.5 m. At 20 h, the typhoon center had left Shantou. However, the water level was still changing, and its distribution pattern was consistent with that in typical weather but nearly 0.2 m higher than in average temperature. In the Chenghai area, near the ocean dam, the outside water level dropped significantly, and the submerged area on the inside remained unchanged. Near the Rongjiang side, the water level outside the dam also dropped, and the inundation inside the dam no longer expanded. The inundation reached its maximum 14 h before, and the deepest accumulated water was 0.38 m. Since then, there was no further disaster expansion.

#### *5.2. Analysis of the Protective Effect of Embankment on Storm Surge*

It can be seen from the storm surge simulation in the 10-year-old situation that the bank can effectively keep out a large amount of seawater. When the typhoon intensity is vigorous, the seawater will cross the bank and enter the urban area. At the same time, after the typhoon transits, the seawater enters the landforms water in the bank instead because of the bank's existence. Figure 9 is a comparison chart of urban inundation degree with or without embankment under the condition of typhoon moving from southeast to northwest during the 10-year return period. The time is the maximum storm surge during the typhoon landing, and the height of the embankment is 1.5 m. It can be seen from the figure that the land in the study area is divided into two parts. One part is adjacent to the South China Sea, and the other is adjacent to the Rongjiang River in Shantou Bay, which can be regarded as the left and right parts of the V-shaped land. When there is no bank, the whole land in the study area will be submerged. The right part of the "V" shape will be submerged seriously, and the accumulated water depth can exceed 0.3 m, while the left part of the "V" shape is higher and far from the sea, and the submerged depth is about 0.1 m. When the bank is set, it can block part of the seawater intrusion, but it will form water near it. The bank can block part of the seawater when the storm surge is most substantial and play a protective role in inland areas.

135

Because the land topography is high, some seawater can flow back into the sea after the storm surge. At the same time, the water increase caused by the storm surge in the ocean will continue to fall back. However, due to the topographic conditions and embankment, part of the seawater will be left on the land. Figure 10 shows the water accumulation on the land when the typhoon is far away from the storm surge. The simulation results show that, in the case of no embankment, the left part of the V-shape is relatively high, and all the backward seawater flows back to the ocean. In contrast, the right part of the V-shape is relatively flat, and the amount of backward seawater is relatively large, so only some seawater can flow back. When there is a bank, the accumulated water on the land cannot enter the ocean due to the block of the bank but accumulates on the inner side of the bank, resulting in higher accumulated water, but less than 0.5 m.

In summary, when a storm surge occurs, the bank can effectively block some seawater from entering the inland and play a better protective role. However, after the storm surge, the seawater entering the land will be blocked by the bank and gathered near the bank, which needs timely dredging.

In order to specifically analyze the blocking effect of the embankment on seawater during the storm surge, the right outside section of the V-shaped embankment was selected, and the changed seawater level and seawater velocity with the connecting distance were analyzed in both conditions. Take the landing time of the typhoon moving southeastnorthwest in the 10-year return period as an example. The comparison of water levels in the sea area outside the embankment is shown in Figure 11. In the figure, a section with a length of 1000 m is selected at the shore bank on the seaside. The red line segment represents the offshore water level when there is a shore bank, and the blue line represents the offshore water level when there is no shore bank. Here, to visually compare the height of the water level, the location of the coastline is taken as the starting surface of the water level, namely 0 m. When there is a bank, the water level near the coastline can reach 1.4 m in this case. The seawater level stays relatively high with the increased offshore distance. However, when the distance is 400–500 m, the water level drops by nearly 0.4 m, which is consistent with the water level without a bank. Without a bank, the water level at the coastline reached 1.03 m under the action of storm surge, which means that the land was submerged by 1.03 m deep seawater at this time. Therefore, it can be said that the embankment significantly impacts the seawater level within 500 m, which is embodied in allowing for the high seawater to gather near the shore.

**Figure 11.** Comparison of water level in the sea area outside the embankment.

The result of comparing the velocity of Shanghai water in the same profile is shown in Figure 12. When there is a bank, because of the barrier of the bank near the coastline, the seawater has no velocity component perpendicular to the shoreline, and the velocity there is small, only 0.21 m/s. With the increased distance from the shore, the velocity of seawater increases in an oscillating way, and the velocity can reach the maximum value of 0.49 m/s at 300 m ashore. After that, the flow rate decreases, reaching the lowest value at 500 m. After 600 m, it becomes stable, and the flow rate is between 0.3 m/s and 0.4 m/s. Under the condition of no bank, the flow rate of seawater in the sea area changes little, and the flow rate on the coastline is faster than that under the condition of no bank. At 600 m, there is a bank gradually approaching the flow rate under the condition of no bank, and at 800 m, the two flow rates are the same, and the changing trend is the same after that.

**Figure 12.** Comparison of velocity in sea area outside embankment.

According to the simulation results, given the typhoon with the intensity of a ten-year return period, the shore dike can effectively block a large amount of seawater when the storm surge occurs. However, due to the function of the shore dike, the offshore water level within 500 m will rise by an additional 0.4 m. When the storm surge is over, the shore dike will cause the seawater entering the land to accumulate water at the shore dike, which needs to be dredged in time. Embankment also has an impact on the velocity of the external sea area within 600 m during a storm surge, resulting in an overall increase of 0.16 m/s in seawater velocity. Generally speaking, the embankment will only work on a specific range of offshore areas.

#### **6. Sediment Deposition in Shantou**

#### *6.1. Simulation Results of Main Siltation Areas*

Figure 13 shows the sediment deposition in the study area over five years. It can be seen that there are four main siltation areas, which are located at Waishahe Estuary, Hanjiang Xixi Estuary, Rongjiang Estuary, and Shantou Port waterway. From the shape of the sedimentary area, the sediments in the three estuaries are all fan-shaped because of the exact sediment deposition mechanism. The Hanxixi Estuary is in a regular fan shape; the Waishahe Estuary and Rongjiang Estuary show an irregular sector, with thicker deposits in the southwest of Waishahe Estuary and vaster deposits in the northeast of Rongjiang Estuary. There are several siltation areas in Shantou Bay, but the siltation degree is lighter than that at the estuary. The siltation in the bay can be divided into two parts: The first is the siltation area near the inland river. The main feature of deposition here is that the deposition range is extensive, but the deposition thickness is small. The second part of the siltation is concentrated on the side near the ocean, showing the trend of cutting off the river. This part of siltation is concentrated and small in area but thick in thickness.

**Figure 13.** Overall siltation in the study area.

Figures 14–17 shows the enlarged views and 5-year changes of siltation thickness of Waishahe Estuary, Hanjiang Xixi Estuary, Rongjiang Estuary, and Shantou Port Channel, respectively, to facilitate the analysis of siltation in the study area. Hanxi River and Waisha River are both branches of the Hanjiang River, so the law of sediment deposition is the same, and the variation curve of sediment thickness in five years is similar.

However, the runoff of Hanjiang Xixi is smaller than that of Waisha River, so the sediment is easier to deposit, and the final deposition thickness is higher than that of the Waisha River Estuary. The central axis of the area where the sediment of Hanxixi Estuary exceeds 1 m is about 1560 m, and the sediment area is about 0.96 km2.

Tides control Rongjiang Estuary and Shantou Bay, and the short-term siltation rates of Rongjiang Estuary and Shantou Bay show periodic changes. A period of 15 days is consistent with the astronomical tide period, indicating that it is highly correlated with the astronomical tide, which is in line with reality. Because of the breakwater in the Rongjiang Estuary, the sediment deposition is different under the ebb and flow tide, forming two settlement centers. The reference point selected by Rongjiang Estuary is located at the sedimentation center outside the sand dike. Different from the Waisha River and Hanjiang Xixi River, the sedimentation rate at the sand dike of Rongjiang Estuary is very uniform and generally increases linearly. Two representative locations in Shantou Port are selected for analysis. The red line segment is the silting area near the sea, and the blue line is the main silting center on the inland river. The sedimentation trend of the two places is the same, but the sedimentation rate near the seaside is faster. The siltation rate in Shantou Bay is decreasing yearly because the siltation in Shantou Bay is more controlled by tides, while the siltation at Rongjiang Estuary and the siltation at the channel reduce the tide carrying sediment into the bay. Therefore, unlike Rongjiang Estuary, the siltation rate in Shantou Bay, which is also controlled by tides, is decreasing yearly.

**Figure 14.** Changes of sedimentation thickness of Waishahe in 5 years.

**Figure 15.** Changes of sedimentation thickness in Hanjiang Xixi Estuary in 5 years.

**Figure 16.** Changes of sedimentation thickness in Rongjiang Estuary in 5 years.

**Figure 17.** Changes of sedimentation thickness in Shantou Bay in 5 years.

#### *6.2. Discuss the Causes of Siltation in Shantou Bay*

Although the siltation range of Shantou Bay is not extensive, the siltation by sea survey will block the river, affect the regular shipping operation, and is not conducive to Shantou's foreign economy and trade. Therefore, it is of far-reaching significance to analyze the causes of siltation in Shantou Bay for the future development of Shantou. Figure 17 shows that the siltation in Shantou Bay is controlled by the tide, which is the decisive factor, but the source of siltation needs further discussion. When rivers flood, the water flow is large, and the sediment content is significant, so the siltation in the study area changes significantly at this time. Therefore, this paper analyzes the influence of rivers on sediment deposition in Shantou Bay from the changes of sediment bodies in Shantou Bay under the historical flood discharge time of each river.

In terms of sediment sources, the sediment sources of Shantou Port mainly include three parts:


Therefore, the historical flood discharge time of the Rongjiang River and Hanjiang Rivers is selected to analyze the siltation areas of Shantou Bay one by one.

For the Rongjiang River, a typical sluice opening and discharging process of the Rongjiang River from 15 August 2016 to 20 August 2016 was selected for analysis. Figure 18 shows the thickness of the siltation body before and after the flood peak here. It can be seen from Figure 18 that after the large-scale flood discharge in the upper reaches of the Rongjiang River, and the maximum siltation area near the inland riverside has increased to some extent. In contrast, the siltation near the seaside has no noticeable change. This shows that the sediment brought by the flood discharge of the Rongjiang River settled before it entered the estuary. It can be judged that Rongjiang's contribution to Shantou Port's siltation will be more negligible if it is not during flood discharge. The siltation contribution of Rongjiang River to Shantou Port mainly occurred during the sluice opening in the upper reaches of Rongjiang River, which mainly acted on the silts near the inland river but made little contribution to the silts near the sea.

**Figure 18.** Changes of silts in Shantou Port before and after flood discharge of Rongjiang River.

The flood discharge of the Hanjiang River and Rongjiang River should be around the 20th of that month to avoid the difference in tidal conditions. Figure 19 shows the comparison of sediment deposition in Shantou Port before and after the flood discharge of Hanjiang River, which is 0:00 on 18 October 2016 and 23:00 on 23 October 2016, respectively. During this period, there is no other flood discharge, and it can be seen that the maximum sedimentation area of the whole sediment body has expanded. Meixi connects Han River and Shantou Port, so Meixi, which directly affects sediment deposition in Shantou Port, is selected for analysis. According to the inspection, the flood discharge of Hanjiang River can reach 670 m3/s, and that of Meixi River is only about 5%, that is, less than 40 m3/s, which is less than 20% of that of Rongjiang River. The discharge of the Meixi River is much smaller than that of the Rongjiang River, but the degree of siltation after flood discharge is greater. Therefore, it can be judged that the Hanjiang River has contributed to the siltation of Shantou Port.

**Figure 19.** Changes of silts in Shantou Port before and after flood discharge of Hanjiang River.

To further discuss the contribution and deposition mechanism of Hanjiang River to Shantou Port, the siltation of Shantou Port before and after the removal of Hanjiang River sediment in one year, and the siltation at Meixi before and after the flood discharge of Hanjiang River in October 2016 were simulated.

Figure 20 to remove the siltation of Shantou Port before and after one year's sediment load of Hanjiang River, the left figure shows the actual siltation situation. The correct figure shows that the sediment load of Hanjiang River is set to 0, and other conditions are consistent to compare the siltation situation of Shantou Port in one year. The sediment deposition in Shantou Port on the right will be reduced to a certain extent, and the sedimentation in the whole area is relatively slow. However, the overall trend has mostly stayed the same. Therefore, it can be concluded that the Hanjiang River contributes to the siltation of Shantou Port, but it is not the decisive factor.

**Figure 20.** Comparison of Hanjiang River's contribution to Shantou Port's siltation. The left picture shows the real state, and the right picture shows that Hanjiang River's silt is set to 0.

Figure 21 shows the change in sediment deposition in the Meixi Estuary before and after the flood discharge. Before the flood discharge in the left picture, it can be seen that there are two prominent silting areas in the Meixi River when it enters the Rongjiang River. In comparison, the silting thickness in the two areas in the right picture has increased after the flood discharge. According to the hydrodynamic conditions of the Hanjiang River and Rongjiang River, the discharge and velocity of the Hanjiang River are much higher. Combined with Figure 21, it can be concluded that the sediment from the Hanjiang River settles after entering the Rongjiang River through the Meixi River, resulting in a deposition at Meixi Estuary. At the same time, the rest continues to be transported downstream. Most of the sediment brought by the single flood discharge of the Hanjiang River settled at the mouth of the Meixi River, and a small part settled in the Rongjiang River.

**Figure 21.** Sediment deposition at Meixi before and after flood discharge.

To sum up, the causes and laws of siltation in Shantou Port are summarized. Sediment sources in Shantou Port include the Rongjiang River, Hanjiang River, and tidal current. Under normal conditions, the Rongjiang River has a small flow rate, slow flow rate, and low sediment load, which does not affect the siltation of Shantou Port. During flood discharge, the sediment of the Rongjiang River mainly settles on the inland river side of Shantou Port. The Hanjiang River has a significant flow rate, a fast flow rate, and a high sediment content, and flows into the Rongjiang River through the Meixi River. However, most of the sediment of the Hanjiang River settles and accumulates at the mouth of the Meixi River after entering the Rongjiang River. Only part of it enters Shantou Port with the Rongjiang River, which contributes to the sedimentation of both parts of Shantou Port, but there are other factors besides this. The decisive factor of sediment deposition in Shantou Port is the sediment entering the port through ocean tides.

#### **7. Conclusions**

The main conclusions of this paper are as follows:

(1) For a typhoon with a 10-year return period, the embankment can effectively block a large amount of seawater in the ocean and prevent seawater from entering the land. At the same time, when the storm surge is over, the seawater entering the ground will not flow back to the ocean smoothly, and water will accumulate at the bank, quickly leading to the aggravation of the post-disaster cleaning task.


**Author Contributions:** Data curation, E.Z. and W.W.; formal analysis, Y.W. (Yuxi Wu) and E.Z.; funding acquisition, Y.W. (Yang Wang) and R.J.; investigation, J.Y., K.Z., H.W. and H.Q.; project administration, Y.W. (Yang Wang) and Z.S.; writing—original draft, Y.W. (Yuxi Wu); writing—review and editing, K.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (Grant No. 52001286; 52101332; 52202427), National Key Research and Development Program of China (Grant No. 2021YFC3101800), Shenzhen Science and Technology Program (Grant No. KCXFZ20211020164015024) and Shenzhen Fundamental Research Program (Grant No. JCYJ20200109110220482).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Wind field data are from the European Centre for Medium-Range Weather Forecasts. Here is the link: https://cds.climate.copernicus.eu (accessed on 11 July 2021).

**Acknowledgments:** European Centre for Medium-Range Weather Forecasts; the China Meteorological Administration.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
