Shape Optimization for the Mitigation of Coastal Erosion via Shallow Water Equations

In October 2003, individuals dug a four-meter breach into the peninsula near the city of Saint-Louis in the Langue de Barbarie, a coastal region in the northern part of Dakar, Senegal. The intention was to counter the flooding that had continually beset the region throughout the years. However, the breach rapidly expanded to 800 meters, separating the southern end of the peninsula from the main country and effectively transforming it into an island. By January 2020, the sea had already claimed more than six kilometers of land; this alteration led to the loss of villages and tourist resorts as well as permanent changes in the peninsula’s flora and fauna. Satellite recordings in Animation 1 capture this process.

While this example is shocking due to its short timescale and resulting severe damage, it highlights a process that will threaten coastal areas all over the world if suitable countermeasures are not taken. The prevention of such disasters amidst continuously rising sea levels is thus limited in neither time nor space, and its importance is only increasing. In order to target the land-displacement process that stems from destructive sea waves, currents, and/or tides, the German Research Foundation currently funds the Shape Optimization Mitigating Coastal Erosion project. This project, which was founded in 2019, supports our collaboration with Diaraf Seck and Mame Gor Ngom from the Cheikh Anta Diop University of Dakar to jointly tackle coastal erosion from a new point of view.

Previous efforts have attempted to mitigate the effects of coastal erosion via groins, breakwaters, and various other structures. Current experiments have utilized the continually increasing computational performance of numerical simulations to model the propagation of waves towards the shore and identify optimal wave-breaking obstacles. Our approach focuses on these kinds of simulations and essentially involves finding an obstacle’s optimal shape for a given wave, relying on shape calculus rather than a finite design space. We hope that this method can yield new and perhaps non-intuitive shapes.

Before targeting the optimization aspect, we must designate a suitable wave description. Hence, we select one of the most widely applied systems of wave equations. We use the Saint-Venant equations—better known as the shallow-water equations (SWE) [2]—to describe the dynamics:

\[\frac{\partial(H)}{\partial(t)} + \frac{\partial Hu}{\partial x}+\frac{\partial Hv}{\partial y} =0\]\[\frac{\partial(Hu)}{\partial(t)} + \frac{\partial Hu^2+ 1/2gH^2}{\partial x}+\frac{\partial Huv}{\partial y} =-gH\frac{\partial z}{\partial x} \tag1\]\[\frac{\partial(Hv)}{\partial(t)} + \frac{\partial Huv}{\partial x}+\frac{\partial Hu^2 + 1/2gH^2}{\partial y} =-gH\frac{\partial z}{\partial y}.\]

The SWE originate from the famous Navier-Stokes equations via depth-integration and are based on the assumption that horizontal length scales are much larger than vertical ones; this notion justifies the terminology of shallowness. It describes the nonlinear time evolution of the water height \(H\) and the weighted velocities \((uH, vH)\) in two dimensions given a sediment height \(z\) (see Figure 1). From a physical point of view, the system maps two conservation laws together: the conservation of mass and momentum, better known as Newton's second law. We often do not need to hold the sediment height as constant in more advanced computations. To exemplify this concept, we can couple the system \((1)\) with an Exner-type equation [4] that describes the conservation of mass for the sediment as

\[\frac{\partial z}{\partial t} = -\frac{1}{1-m} \nabla \cdot \mathbf{Q} \tag2\]

for porosity \(m\) and sediment flux \(\mathbf{Q}\). We can then formulate a constrained optimization problem in the simplest form, where flux matrix \(F\) and sources \(S\) are defined in accordance with \((1)\)-\((2)\):

\[\min J(\Omega, \mathbf{U}) = \int_0^T \int_{\Gamma_\textrm{Shore}} \frac{1}{2}[H(t,x) - \tilde{H}(t,x)]^2 \textrm{d}s \, \textrm{d}t+v_1 \int_D 1 \textrm{d}x \tag3\]\[\textrm{s.t.}  \enspace \enspace \partial_t\mathbf{U}+\nabla \cdot F(\mathbf{U})=S(\mathbf{U}) \enspace \textrm{on} \enspace \Omega. \]

We can interpret \((3)\) as the search for a shape of an obstacle \(D\) that minimizes the squared difference of \(H\) to a target value \(\tilde{H}\) at the shore \(\Gamma_\textrm{Shore}\) for solution \(U=(H, uH, vH, z)\) at sea level \(\Omega\). Some suitable sub-goals can accompany this objective; here we display a volume penalty that is controlled by \(v_1\). In addition, recent work extends this penalty by a perimeter regularization and a thickness control [6].

Identifying the extreme values of a function generally requires knowledge of derivatives. When looking at the Lagrangian of \((3)\)—as is common in equality-constrained optimization—we are interested in a shape analogue to the classical derivative, wherein we build upon existing ideas [3, 8]. We obtain such a derivative by evaluating the Lagrangian in the saddle point; doing so allows us to determine the solution to the coupled SWE \((1)\)-\((2)\) as well as the so-called adjoint equations. The solution to both problems is required when computing the shape derivative, which we ultimately use to deform the domain \(\Omega\).

Since analytical solutions to \((3)\) are unavailable, we solve this problem numerically with a gradient descent algorithm. The development of numerical solutions to the SWE is an active field, and numerous studies have investigated classical schemes like finite difference, finite element, and finite volume methods. However, all approaches generally calculate a discretized solution on a partition—commonly called a mesh—of the domain \(\Omega\). We instead base our computations on a discontinuous Galerkin scheme that inherits advantages from both finite volume and finite element methods. 

Revisiting the Langue de Barbarie once more, we hence need to construct a mesh that represents the domain, wherein we calculate a solution to the partial differential equation at each mesh point. We utilized the free Global Self-consistent, Hierarchical, High-resolution Geography Database to obtain the mesh points and created the mesh in Figure 2 with the finite element mesh generator Gmsh [5]. In addition, we consulted the General Bathymetric Chart of the Oceans database to map the real sediment height to the mesh points via a nearest-neighbor classification.

In the setting of \((3)\), we used FEniCS to simulate a tsunami-like water evaluation for an open-sea boundary at the half circle and no-normal flow boundary conditions at an island, shore, and obstacle [1]. Figure 3 illustrates the propagating wave at different time points.

We target a rest height at \(\Gamma_\textrm{Shore}\), as pictured in Figure 2. Initially, we placed a circular obstacle on the mesh. While fixing the remaining boundaries, we optimized the obstacle's shape using the aforementioned algorithm. We were able to achieve convergence to a new shape and drastically minimize the objective. The results are evident in Figure 4.

Because the obstacles that we obtain are often too large for practical implementations—e.g., the initial diameter of the circular obstacle in the calculated example is roughly 500 meters—we acknowledge that this work can only serve as an initial feasibility study. We therefore hope to increase the realism of the simulations in the near future. Nevertheless, we have already extended the described approach to arbitrary coastal sections and multiple obstacles — such as various initial and boundary conditions for different objectives. We have also computed results on spherical world meshes as solutions on immersed manifolds, which pave the way for global simulations.

Two associated working papers are currently available and are updated on a regular basis [6, 7]. 

Luka Schlegel presented this research during a minisymposium at the 2021 SIAM Conference on Control and its Applications, which took place virtually in July 2021 in conjunction with the 2021 SIAM Annual Meeting.

Acknowledgments: This work was supported by the German Research Foundation within the priority programme SPP 1962, "Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and Hierarchical Optimization." The authors would like to thank Diaraf Seck and Mame Gor Ngom of the Université Cheikh Anta Diop of Dakar for helpful and interesting discussions within the Shape Optimization Mitigating Coastal Erosion project.

References

[1] Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., …, Wells, G.N. (2015). The FEniCS project version 1.5. Archive Num. Software, 3(100).
[2] Barré de Saint-Venant, A.-J.-C. (1871). Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et è l'introduction des marées dans leur lit. C. R. Acad Sci Paris. 
[3] Delfour, M.C., & Zolésio, J.-P. (2011). Shapes and geometries: Metrics, analysis, differential calculus, and optimization. Philadelphia, PA: Society for Industrial and Applied Mathematics. 
[4] Exner, F.M. (1925). Über die Wechselwirkung zwischen Wasser und Geschiebe in Flüssen. Sitzungsberichte. Akademie der Wissenschaften in Wien.
[5] Geuzaine, C., & Remacle, J.-F. (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng., 79(11), 1309-1331. 
[6] Schlegel, L., & Schulz, V. (2021). Shape optimization for the mitigation of coastal erosion via shallow water equations. Preprint, arXiv:2107.09464.
[7] Schlegel, L., & Volker, S. (2021). Shape optimization for the mitigation of coastal erosion via the helmholtz equation. Preprint, arXiv:2107.10038.
[8] Sokolowski, J., & Zolesio, J.-P. (1992). Introduction to shape optimization: Shape sensitivity analysis. In Springer series in computational mathematics. New York, NY: Springer-Verlag.