1 Introduction

Coral reefs are representative components of marine geomorphology, are widely distributed in the deep sea and shallow waters, and especially exist in warm waters on both sides of the equator (Putnam et al., 2017). In particular, coral reef topography is widely distributed in Taiwan Island, the South China Sea, and Australian waters (Pirazzoli et al., 1993; Harris et al., 2004; Beaman et al., 2008; Woodroffe and Webster, 2014). In recent years, with the rapid expansion of marine resources and the ongoing construction of civil structures and military facilities on coral reefs, there is an urgent need to carefully study the stability of coral reef slopes. The regularity of horizontal distribution of coral reefs tends to form steep seaward slopes, hard reef flats and coral sand lagoons covered with hard reef shells depending on the hydrodynamic environment and geological conditions, as well as the long-term formation process of coral reefs in the South China Sea (Rogers et al., 2017). In this special marine environment and geological condition, a typical two-layer profile appears on the vertical cross-section, with incompact coral sand to a depth of 10 - 20 m consisting of coral debris, gravel, and coarse sand; and a deep layer of hard coral reef with good integrity (Gong et al., 2012; Kordi and O'Leary, 2016). A typical profile of the geomorphic features of coral reef in the South China Sea is shown in Fig.1 (Zhang et al., 2017), wherein the coral reef is covered with a hard shell at the surface. As shown in Fig.1, the geological structure of the coral reef consists of the following parts: 1) the upper seaward slope; 2) the outer reef flat; 3) the reef margin; 4) the inner reef flat; 5) the lagoon slope; 6) the lagoon basin. The first four parts (1 - 4) are strongly subjected to hydrodynamic action.

Such projects as the construction of ports and military structures, along with shipping and tourism activities, tend to destroy the hard shell of coral reefs, which may significantly reduce their general stability; in turn, this can reduce the stability of structures established on the reefs (Chen et al., 2017). Therefore, the integrity of the reef shell has a great influence on ensuring the stability of coral reefs. In order to investigate theoretically the influence of horizontal wave loads on the stability of coral reefs, and in light of energy conservation considerations, we propose an energy functional for coral reef stability resulting from the combined vertical self-gravity and horizontal wave loads. Then, we obtain a variational solution for the stability of coral reefs covered with a hard shell under horizontal wave loading, thereby demonstrating the influence of the horizontal wave loads and the original reef shell on the reef slope stability. We simplified the geomorphic features of the coral reef in order to reduce the difficulty of theoretical derivation and mathematical calculation, as shown in Fig.2. Table 1 presents data on the hydrodynamic and geological conditions related to the coral reef, obtained from the actual monitoring of the sea area and the coral reefs.

According to a previous study, the variational limiting equilibrium (LE) solutions are identical to the upper plastic limit analysis solutions; thus, the top hard shell can obviously improve reef stability; however, this study ignored the horizontal action due to the wave loads (Zhang et al., 2018). The main advantage of this variational approach is that it is free from any artificial kinematical or static admissible assumptions, which is the case for many existing methods mentioned in the literature. In order to assess the stability of coral reefs under horizontal wave loading, the concept of safety factor is adopted in the current paper, which is usually defined as the ratio of the total shear strength available on the slip surface to the total shear strength required for equilibriums.

2 Variational Analyses

Here, we theoretically derive the stability factor of the coral reef using the variational LE in order to determine the reef stability covered with hard reef shells and undergoing the horizontal wave loads. The analysis process is as follows.

2.1 Mathematical Formulation of the Basic Problem

A typical reef slope i inclined in the horizontal direction is shown in Fig.3. The coral reef profile is identified by the unit weight γ, cohesion c, and internal friction angle φ. Here, h is an assumed height of coral reef described in this paper. The horizontal wave load p(x), the vertical load on top of the reef Q(x), and the tension crack ζ on the reef roof are introduced, assuming that the collapse of the coral reef is due to the combined vertical self-gravity load and horizontal wave loads. The LE method is based on the following equations:

1) Satisfaction of the plastic yielding criterion τ = f(σ) along the shear slip surface y(x). We adopt the Mohr-Coulomb's failure condition as:

$ \tau = c + \sigma \psi, $

(1)

where τ = τ(x) and σ = σ(x) are the distributions of the tangential and normal stress along the assumed slip surface y(x), respectively, and D_s refers to the thickness of the surface reef shell. To facilitate easy calculation in this paper, we set ψ = tan(φ).

2) Satisfaction of the force and moment equilibrium equations for the sliding mass bounded by slope surface $ \bar y(x) $ and shear slip surface y(x).

3) Here, the horizontal wave load p(x) is applied on the coral reef, and the wave shape adopted before wave breakage is calculated as in Fig.4. Considering the complexity of wave force calculation theory on slope, and in order to simplify the theoretical derivation of horizontal wave pressure, it is assumed that the reef shell is inclined at a fixed angle i to the horizontal plane; besides, the surface roughness of the reef shell and the wave reflection effect are ignored.

The wave pressure p (kPa) acting on the coral reef is calculated as follows:

$ p = 10{k_1}{k_2}{\bar p_2}\gamma H, $

(2a)

where k₁ = 1.35 and

$ {k_2} = 0.85 + \frac{{4.8H}}{L} + \cos i(0.028 - \frac{{1.15H}}{L}), $

and they are the coefficients; ${\bar p_2} = 1.7$ is the relative maxi-mum wave pressure ${\bar p_2} = 1.7$; the parameters H and L are the wave height and wave length, respectively; and i is the slope inclination of the coral reef.

The vertical coordinates z₂ and z₃ of point 2 and point 3 are respectively determined by the following equations:

$ {z_2} = A + \frac{1}{{\cos i}}\left[ {1 - \sqrt {2{{\cos }^2}i + 1} (A + B)} \right], $

(2b)

$ {z_3} = H\left[ {\frac{{\sqrt {\frac{{\text{π }}}{{2i}}} + {\text{π }}\frac{H}{L}\coth 2{\text{π }}}}{L}d} \right], $

(2c)

where A and B are the formula parameters, which can be found in the Soviet architecture standard.

Based on the coral reef mass being in a state of LE, as shown in Fig.3, the equilibrium equations of horizontal and vertical forces and the moment are respectively satisfied as follows (Zhang et al., 2008):

$ H = \int_s {\left({\tau \cos \alpha - \sigma \sin \alpha } \right){\text{d}}s} + \int_{{x_0}}^{{x_2}} {p{\text{d}}x} = 0, $

(3a)

$ V = \int_s {\left({\tau \sin \alpha + \sigma \cos \alpha } \right){\text{d}}s} -\\ \int_{{x_0}}^{{x_2}} {\gamma \left({\bar y - y} \right){\text{d}}x} - \int_{{x_0}}^{{x_2}} {Q(x){\text{d}}x} = 0, $

(3b)

$ M = \int_s {\left[ {\left({\tau \sin \alpha + \sigma \cos \alpha } \right)x - \left({\tau \cos \alpha - \sigma \sin \alpha } \right)y} \right]{\text{d}}s} -\\ \int_{{x_0}}^{{x_2}} {\left[ {p\bar y + \gamma \left({\bar y - y} \right)x} \right]{\text{d}}x - \int_{{x_0}}^{{x_2}} {Q(x)x{\text{d}}x} }, $

(3c)

where s is the arc length along the shear sliding surface y(x), and α is the slope inclination of y(x). Eqs. (3a) and (3b) are the horizontal and vertical force equilibrium equations in the two-dimensional (2D) plane strain state, respectively, and Eq. (3c) is the moment equilibrium equation of the rigid-plastic sliding coral reef. We then introduce Mohr-Coulomb's failure criteria into the above equations and combine the geometrical relations along the continuous boundaries shown in Fig.3, cosα = dx/ds; sinα = y'dx/ds. Under these simplifications, we can obtain the following equilibrium equations which can indicate the overall stability equations of this layered slope:

$ H = \sum\limits_{i = 0}^2 {{h_i}} = \sum\limits_{i = 0}^2 {\int_{{x_i}}^{{x_{i + 1}}} {\left[ {{c_i} + {\sigma _i}\left({{\psi _i} - y'} \right) + p} \right]{\text{d}}x} } = 0, $

(4a)

$ V = \sum\limits_{i = 0}^2 {{v_i}} = \sum\limits_{i = 0}^2 {\int_{{x_i}}^{{x_{i + 1}}} {\left[ {{\sigma _i}\left({{\psi _i}y' + 1} \right) + {c_i}y' - {\gamma _i}\left({\bar y - y} \right) - Q} \right]} } {\text{d}}x = 0, $

(4b)

$ M = \sum\limits_{i = 0}^2 {{m_i}} = \sum\limits_{i = 0}^2 {\int_{{x_i}}^{{x_{i + 1}}} {\left[ {{c_i}{x_i}y' + {\sigma _i}{x_i}\left({{\psi _i}y' + 1} \right) - {c_i}y - {\sigma _i}y\left({{\psi _i} - y'} \right) - {\gamma _i}\left({\bar y - y} \right){x_i} - Qx - p\bar y} \right]} {\text{d}}x = 0} . $

(4c)

2.2 Safety Factor and Critical Height of Reef Slope

To quantify the margin of reef safety relative to an LE state under combined horizontal and vertical loads, we replace the real shear strength parameters of coral reefs with an artificially reduced factor F. This is also known as safety factor and is expressed as:

$ \bar c = c/F, \bar \psi = \psi /F . $

(5)

The critical height h_cr of the coral reef under combined vertical and horizontal wave loads, for which we can obtain the reef stability when it reaches a state of LE, is decided by the normal stress functions and kinematical shape σ(x) and y(x) along the shear sliding surface, respectively. Thus, the reduced factor F and slope height h can be described as a variational functional of the mentioned shape and stress functions. Consequently, the safety factor F_s and the critical height h_cr are respectively considered the extreme values, as follows:

$ {F_s} = \mathop {\min }\limits_{\left({y, \sigma } \right)} F\left\{ {y\left(x \right);\sigma \left(x \right)} \right\}, $

(6a)

$ {h_{cr}} = \mathop {\max }\limits_{(y, \sigma)} h\{ y(x);\sigma (x)\} = h[{y_{cr}}(x), {\sigma _{cr}}(x)], $

(6b)

where y_cr(x) and σ_cr(x) are the shear sliding surface and normal stress distribution along the optimal failure slip surface, respectively.

In the process of solving the functional extremum shown in Eq. (6), we set the safety factor F_s to 1 and use the critical height h_cr to estimate the stability of the coral reef. Thus, the LE issue can be stated as solving two sets of mathematical functions y(x) and σ(x), which can obtain the extreme value h_cr of the functional and simultaneously satisfy the coral reef equilibrium equations (Eq. (4)) mentioned above.

We can transform the basic equilibrium issues into the standard isoperimetric problem of the variational method using the derivation process previously introduced by Baker(1981, 2003). We select the force equilibrium equation in the vertical direction as the integral object, using the horizontal force and moment equilibrium equation as the integral constraints. The parameters λ₁ and λ₂ are Lagrange's undetermined multipliers, respectively, which are constants in the process of variational derivation. Then, by introducing the two integral constraints (Eqs. (4a) and (4c), respectively), we can obtain the following expression g[y(x), σ(x)] given by

$ g = \sum\limits_{i = 0}^2 {{g_i}} = \sum\limits_{i = 0}^2 {\int_{{x_i}}^{{x_{i + 1}}} {{v_i} + {\lambda _1}{h_i} + {\lambda _2}{m_i}} } = \sum\limits_{i = 0}^2 {\int_{{x_i}}^{{x_{i + 1}}} {\left\{ {{\sigma _i}\left({{\psi _i}y' + 1} \right) + {c_i}y' - {\gamma _i}\left({\bar y - y} \right) - Q + {\lambda _1}\left[ {{c_i} + {\sigma _i}\left({{\psi _i} - y'} \right) + p} \right]} \right.} } $

$ \left. { + {\lambda _2}\left[ {{c_i}{x_i}y' + {\sigma _i}{x_i}\left({{\psi _i} + 1} \right) - {c_i}y - {\sigma _i}y\left({{\psi _i} - y'} \right) - {\gamma _i}\left({\bar y - y} \right){x_i} - Q - p\bar y} \right]} \right\}{\text{d}}x . $

(7a)

$ \mathop {\max }\limits_{(y, \sigma)} g = 0 . $

(7b)

We have to obtain the optimal sliding surface and its corresponding normal stress distribution in order to solve the extremum of the above auxiliary functional shown in Eq. (7b).

2.3 Functional Variational Solution

In order to solve the safety factor F_s and critical height h_cr of the coral reef, the mentioned functional must satisfy the Euler partial differential equations, i.e., the necessary conditions for the existence of an extreme, namely,

$ \frac{{\partial {g_i}}}{{\partial \sigma }} - \frac{{\text{d}}}{{{\text{d}}x}}\left[ {\frac{{\partial {g_i}}}{{\partial \sigma '}}} \right] = 0, $

(8a)

$ \frac{{\partial {g_i}}}{{\partial y}} - \frac{{\text{d}}}{{{\text{d}}x}}\left[ {\frac{{\partial {g_i}}}{{\partial y'}}} \right] = 0 . $

(8b)

According to Baker and Garber(1977, 1978), we can derive the family of the potential shear sliding surface by solving the first Euler partial differential equation, whereas the second Euler equation could yield the normal stress acting along the shear sliding surface. In this paper, the shear slip surface equation and the normal stress equation are derived easily by solving the Euler partial differential equations (Eqs. (8a) and (8b)), which could be expressed as Eqs. (9a) and (9b). The derivation process is consistent with that followed by Baker. The final theoretical formulas of the slip surface and its corresponding normal stress are as follows:

$ \rho = \left\{ {\begin{array}{*{20}{c}} {{\rho _o}\exp \left({ - \psi \beta } \right), {\text{ }}\psi \ne 0} \\ {{I_1}, {\text{ }}\psi = 0} \end{array}} \right., $

(9a)

$ \sigma = \left\{ {\begin{array}{*{20}{c}} {\Omega \cdot \exp \left({2\psi \beta } \right) + \frac{{\gamma {\rho _o}}}{{1 + 9{\psi ^2}}}\exp \left({ - \psi \beta } \right)\left({\cos \beta + 3\psi \sin \beta } \right) - \frac{c}{\psi }, {\text{ }}\psi \ne 0} \\ {2c\beta + \gamma \rho \cos \beta + {I_2}, {\text{ }}\psi = 0} \end{array}} \right., $

(9b)

where ρ_o, Ω, I₁, and I₂ are the integration constants. Note that Eqs. (9a) and (9b) are formulations of the family of the potential sliding surface and its corresponding normal stress in a polar coordinate system, respectively. We utilized the following coordinate transformation, which can convert the coordinate values from polar to rectangular coordinates:

$ x = - \frac{1}{{{\lambda _2}}} + \rho \sin \beta = {x_c} + \rho \sin \beta, $

(10a)

$ y = \frac{{{\lambda _1}}}{{{\lambda _2}}} - \rho \cos \beta = {y_c} - \rho \cos \beta, $

(10b)

where (ρ, β) is the polar coordinate system with a center (x_c, y_c), as shown in Fig.3.

Next, we can integrate the force and moment equilibrium equation explicitly once the form of the potential sliding surface and its corresponding normal stress distribution are obtained. The final form of this equilibrium equation is described as

$ H = V = M = 0 . $

(11)

It is important to note that we must still consider the constraints listed below in order to obtain a specific theoretical solution of the above Euler partial differential equation.

1) Geometrical boundary conditions

$ y\left({x = {x_0}} \right) = \bar y\left({x = {x_0}} \right) = 0, $

(12)

$ y\left({x = {x_2}} \right) = {h_{cr}} - \zeta . $

(13)

Note that the two endpoints, 0 and 2, of the sliding surface are not fixed, i.e., they are separately variable along the slope surface. The tensile crack depth ζ can be obtained using the tensile strength of the coral reef. The longitudinal coordinates are already given.

2) Transversality condition among layers

In the process of searching for the optimal sliding surface, the coordinates of the two endpoints, 0 and 2, cannot ascertain which are related to the loading and the geological properties, so it is necessary to apply the transversality conditions among the layers, whose general form is given below:

$ {\left. {\left[ {g + \left({\Theta ' - y'} \right)\frac{{\partial g}}{{\partial y'}}} \right]} \right|_{x = {x_i}}} = 0, $

(14)

where x_i may be either x₀, x₁, and x₂. The Θ is a particular curve passing through the two endpoints. By substituting relevant expressions into Eq. (14), utilizing the coordinate transformations in Eq. (10), and then making a detailed calculation, the following expression is obtained:

$ {\left. {\left[ {\sigma \left(\beta \right)\left({\sin \beta + \psi \cos \beta } \right) + c\cos \beta - \gamma \left({\bar y - y} \right)\sin \beta } \right]} \right|_{{\beta _i}}} = 0 . $

(15)

By applying the above Eq. (15) to the endpoint 0, we obtain the following expression:

$ {\sigma _1} = \frac{{ - {c_s}}}{{\tan {\beta _1} + \psi }}, $

(16)

where the parameter c_s indicates the cohesion of the coral reef shell.

Simultaneously, we apply the equation to endpoint 2 and consider the existence of tensile crack ζ, that is, $ (\bar y - y) $= ζ. On the one hand, if the tensile crack ζ on the coral reef surface is well-known, then we can obtain the following expression:

$ {\sigma _2} = \frac{{ - {c_s} + \gamma \zeta \tan {\beta _2}}}{{\tan {\beta _2} + \psi }} . $

(17a)

On the other hand, if the ultimate tensile strength T of the hard reef roof is known, we can obtain the depth of tension crack ζ, as follows:

$ \varsigma = \frac{{\tan {\beta _2} + \psi }}{{\gamma \tan {\beta _2}}}T + \frac{c}{{\gamma \tan {\beta _2}}} . $

(17b)

3) Stress boundary condition

Let us substitute the stress boundary condition σ₂ = σ(β = β₂) into Eq. (9). In doing so, we can solve the integral constants Ω and I₂, after which the following expressions are obtained as shown in Eq. (18):

$ \Omega = \left[ {{\sigma _2} - \frac{{\gamma {\rho _o}}}{{1 + 9{\psi ^2}}}\exp (- \psi {\beta _2})\left({\cos {\beta _2} + 3\psi \sin {\beta _2}} \right) + \frac{c}{\psi }} \right] \cdot \exp \left({ - 2\psi {\beta _2}} \right), {\text{ }}\psi \ne 0, $

(18a)

$ {I_2} = - \frac{c}{{\tan {\beta _2}}} - 2c{\beta _2} - \gamma \rho \cos {\beta _2}, {\text{ }}\psi = 0 . $

(18b)

Hence, we can statically determine the basic variational problem by simultaneously solving the equilibrium equations as well as the geometrical boundary, stress boundary, and transversality conditions. The following two calculation results, which consider different slope inclinations and reef tensile strengths shown in Fig.5, give the typical shear sliding surfaces and their corresponding normal stress distributions with slope angles π/4 and π/2, respectively.

2.4 Numerical Verification of Theoretical Solutions

As shown in Fig.5, the general finite element software Plaxis is used to simulate the engineering example. This is done in order to verify the rationality of the variational results derived theoretically in this paper and to check the shape of the sliding surface and the normal stress distribution law along the sliding surface. The numerical calculation results are shown below.

The shear sliding surfaces shown in Figs.6(a) and (c) indicate that the hard surface shell constrains the shallow sliding of the coral reef and increases the general stability of the coral reef. The safety factor N_s of the coral reef with four differential potential critical sliding surfaces is shown in Fig.6(a), which indicates that the sliding surface gradually transits from shallow failure to deep failure, increases from 16.6 to 28.6. In addition, Fig.6(c) shows that the safety factor increases from 4.16 to 10.3. Furthermore, the failure modes shown in Figs.6(b) and (d) indicate that the appearance of tensile cracks at the top surface of the coral reef accelerates the appearance of shear sliding surface and reduces the general stability of the coral reef. Furthermore, the failure modes of the numerical calculation are consistent with those of the variational solution derived in this paper. Therefore, protecting the integrity of the coral reef shell is very important in ensuring the overall stability of the coral reef.

3 Comparison and Discussion of Results 3.1 Variational Solution of the Stability Factor N_s

Actually, it is difficult to obtain the analytical results of the nonlinear simultaneous equations derived theoretically in this paper. For this reason, a trial-and-error process is utilized to generate the numerical solutions. To facilitate easier calculation, we express all the quantities in their non-dimensional forms. For example, the length-type parameters take the expression of l = γL/c, and the stress-type parameters are expressed as s = S/c. Furthermore, we use a non-dimensional parameter, which has been defined as N_s = γh_cr/c by Taylor (1948). By doing so, we obtained some theoretical results indicating the critical safety factors for various conditions with the different slope inclinations as well as the internal friction angle of the coral reef, as illustrated in Fig.7. By comparison, we find that the safety factor derived using the variational LE procedure is similar to that obtained by the upper bound analysis method.

Fig.7 shows that the safety factor N_s is a function of φ and i (internal friction angle and reef inclination, respectively) for the case of no reef shell covered. Comparing the results with the upper bound analysis results published by Chen (1975), we can expect the variational approach to obtain the exact theoretical values of N_s, and to derive a set of theoretical calculation fourmulas.

3.2 Safety Factor with Hard Reef Shell

Some investigators have concluded that the influence of shell thickness D_s on the coral reef stability (i.e., to cover the surface of a coral reef) is not negligible. Fig.8 presents N_s as a function of internal friction angle φ and reef inclination i for the case of a coral reef with a hard shell (D_s). On the one hand, the dotted lines in Fig.8 represent the reef safety factor without the reef shell being covered. On the other hand, the straight lines mean that with a hard reef shell, the existence of a continuous covered reef shell significantly increases the stability of the coral reef.

3.3 Safety Factor with Tensile Crack

Fig.9 shows the effect of surface tensile strength T on the safety factor N_s with the emergence of the tensile crack ζ. We can see that there exists a linear growth relationship between the safety factor N_s and the tensile strength T. Furthermore, the figure shows that the safety factor N_s of the coral reef is approximately linearly related to the surface tensile strength T, and the line slopes are basically the same under different reef inclination angles i.

3.4 Safety Factor Due to Horizontal Wave Loads

Two sets of typical wave data (Figs.10 and 11), which include water depth in front of the coral reef 40 m, wave period 4 s, wave length 24.98 m, and wave heights 1 and 5 m, respectively. The wave loads p(x) acting on the reef are solved, and the safety factor N_s of coral reefs is analyzed under wave loads and self-gravity. Based on the variational calculation results, we determine that the safety factor N_s of the coral reef decreases with the increase of wave load p.

4 Conclusions

In this paper, the variational LE procedure is applied to obtain the safety factor N_s of the coral reefs covered with hard reef shells and to withstand wave loads in the horizontal direction. Some meaningful conclusions are drawn as follows:

1) The analysis results show that the continuous hard shell covering the coral reef enhances the stability of the coral reef. Furthermore, there is a nonlinear relationship between the safety factor N_s and the tensile strength T and the thickness D_s of the surface reef shell.

2) The horizontal wave loads weaken the safety factor N_s of the coral reef, and the existence of tensile crack ζ further reduces the safety factor N_s of the coral reef.

Acknowledgements

This paper is supported by the Project of National Science and Technology Ministry (No. 2014BAB16B03), and the National Natural Science Foundation of China (No. 51679224).