Three-Dimensional Seismic Response Analysis of Zhubi Reef in the South China Sea

Reef islands are valuable terrestrial resources for ocean exploitation and utilization but face the serious threat of earthquake disasters. *is study takes Zhubi Reef in the South China Sea as the research object and establishes a three-dimensional (3D) reefseawater model to investigate the seismic response of Zhubi Reef.*e 3D local topography and the effect of fluid-solid interaction are comprehensively considered. *e artificial boundaries of fluid and solid domains are adopted to simulate the wave radiation and absorption of the semi-infinite seabed and the infinite seawater layer. *e boundary substructure method is used to input the seismic waves into the numerical reef models. *e results indicate that the seismic responses, including the peak values of the ground motions and the acceleration response spectra, are significantly amplified over Zhubi Reef. *e acceleration response spectra on the reef flat shift closer towards the low-period direction compared with those of the seismic waves input from the subsea bedrock. In addition, the 3D topographic effect on the seismic response of Zhubi Reef is studied through a comparative analysis between two-dimensional (2D) and 3D reef models.*e distribution laws of the peak response in the middle region of the reef flat calculated by the 2D and 3Dmodels agree well with each other, and the differences are obvious in the edge areas where the local 3D topographies change drastically.


Introduction
With technological progress and economic development, the exploitation and utilization of the ocean have entered a new stage. Reef islands are valuable terrestrial resources in the sea. In recent years, large-scale artificial islands have been constructed based on islands and reefs in the South China Sea, on which a series of important infrastructures, such as airport runways, port terminals, and large lighthouses, have been built. ese constructions provide important footholds and supply depots for activities such as fishery production, shipping transportation, scientific research, and meteorological observation in the sea area. However, the island and reefs in the South China Sea are located on the west side of the circum-Pacific seismic belt, which is a vulnerable zone facing the serious risk of earthquake disasters [1]. Although researchers have paid more attention to the engineering problems of reef sites, research on the dynamic response of reef islands and the key ground motion parameters of reef flats is still insufficient to complete the seismic design of the artificial islands and the reef constructions.
Previous studies on site seismic responses were mostly aimed at land-based engineering sites [2][3][4][5][6][7]. Abundant research resources in observation data, analysis methods, and response mechanisms and regularities have been accumulated for this research field. Reasonably improving the relevant land-based site seismic analysis methods and applying them to the seismic analysis of islands and reefs have become the main research focus for current seismic resistance studies of reef engineering. Based on this idea, Hu et al. [8] studied the amplification effect of coral reefs on pulse seismic waves through a one-dimensional seismic response analysis method of layered soil sites. en, they established a 2D reef model with a height of 40 m, using equivalent nodal forces to simulate seawater but ignoring the dynamic interaction between fluid and solid, to study the seismic response of reef islands [9]. Chen et al. [10] considered the nonlinear dynamic characteristics of coral sand and the influence of the infinite foundation, established a 2D seismic analysis model of coral reefs with a height of 20 m, and analyzed the spatial distributions of the ground motion field and the response spectrum on reef flats. However, the effect of fluid-solid interactions is not considered in this study.
In fact, there are several obvious differences in the seismic analysis models and response mechanisms between reef engineering sites and traditional land-based sites, which can be summarized as follows. (1) e topographic and geographic characteristics of coral reefs are unique [11,12]. Reef islands are deposited by the remains of corals from the deep water environment of the South China Sea at a depth of 1400-2000 m throughout an extremely long geological period. Artificial reefs are constructed with coral sand. e topography and geomorphology of reef islands as well as the physical and mechanical properties of coral reefs and sand have an important influence on the seismic response law and the distribution of ground motion fields of reef engineering sites. (2) Under the dynamic coupling effect of infinite seawater and reef islands, due to the large scale of the reef body, compression waves in seawater will have an important effect on the dynamic response of the reef. erefore, the compressibility and wave radiation effect of seawater should be considered in the analysis model. In addition, the interaction between the reef islands and the semiinfinite seabed is also an important issue when calculating the seismic response of a reef model. e key techniques involved in this problem include the simulation of a semi-infinite medium and the input of seismic waves. In the previous studies completed by our research team, we comprehensively considered the abovementioned key factors and established a seismic analysis model of reef islands [13]. e artificial boundaries of fluid [14] and solid [15] media are applied on the truncated boundary of the near-field model to simulate the wave radiation effect, and the seismic wave input method based on the substructure of artificial boundaries (the boundary substructure method, BSM) [16] is adopted to realize the incidences of seismic P and SV waves at different angles in the reef models. en, a typical reef island in the South China Sea, namely, Yongshu Reef, is taken as a research object, and its ground motion field distribution and seismic response characteristics are studied through 2D and 3D analysis [17,18]. e reef islands in the South China Sea can be mainly divided into two types according to their morphologies. e first type, i.e., the Yongshu Reef, presents long and narrow shapes, and the other is the atolls represented by Zhubi Reef. Compared with Yongshu Reef, Zhubi Reef has a larger reef flat and thus a greater potential for engineering development. On account of its enormous geometric scale, however, there are great challenges in the 3D seismic response analysis of Zhubi Reef. is study takes Zhubi Reef as the research object, comprehensively considers the 3D local topographic characteristics, the effect of fluid-solid interaction, the wave radiation effect of the semi-infinite seabed and the infinite seawater layer, and the seismic wave input method, and establishes a 3D reef-seawater system model to investigate the distribution of the seismic wave field on Zhubi Reef. In addition, the 3D topographic effect on the seismic response of Zhubi Reef is studied through a comparison analysis between the 2D and 3D reef models.

Numerical Models and Seismic Analysis
Methods of Reef Sites

Numerical Model and Model
Parameters. e local map of Zhubi Reef in the South China Sea is shown in Figure 1.
e reef flat presents a pear shape. e length and width of the reef flat are 5.8 km and 3.4 km, respectively [19]. Based on the geometrical characteristics of the reef flat and the collected material parameters of Zhubi Reef, we establish a 3D calculation model of the reef site, as shown in Figure 2.
e height of the reef body in this model is 800 m, and the slope gradients are basically distributed between 15°and 30°. e depth of the seabed intercepted from the semi-infinite domain is 400 m, and the minimum width of the surrounding seawater, namely, the distance from the foot of the reef to the truncation boundary of the sea region, is 1000 m. e reef body is divided into three layers along the vertical direction, which are denoted as layers 1, 2, and 3. e material parameters of the reef-seawater model were determined according to Shan et al. [20] and Sun and Lu [21] and are presented in Table 1. Considering that the depth of the lagoon on Zhubi Reef is relatively shallow (less than 20 m) [19], the effect of the lagoon on the dynamic response of Zhubi Reef is ignored in this study. Six typical slices are extracted from the 3D reef model along the northwest-southeast direction, whose lengths of reef flats are 2020 m, 2935 m, 3140 m, 2930 m, 1500 m, and 1010 m. To evaluate the 3D topographic effect on the seismic response, six corresponding 2D reef models with identical topographic and material parameters as those of the 3D model are established, as shown in Figure 2. e general finite element program ANSYS is adopted to establish the reef-seawater models, as shown in Figure 3. e compressible acoustic fluid element Fluid80 and the solid hexahedral element Solid45 are applied to discretize the seawater region and the reef-seabed regions, respectively. e reef-seawater interactions are simulated by coupling the normal degrees of freedom of the fluid and solid nodes in the interfaces [13].

Fluid and Solid Artificial Boundaries.
e stress-type fluid artificial boundary [14], which converts the partial differential equation of the wave motions in fluid into spatially decoupled equivalent damper-mass systems, is adopted to absorb the outgoing scattered waves generated in the seawater region, as shown in Figure 3. Under 2D and 3D conditions, the mass and damping coefficients of the fluid artificial boundary are given as follows: where ρ F is the mass density of fluid; R F is the distance from the wave source to the cutoff boundary in the seawater domain; c F is the sound velocity of fluid; and Σ i A i stands for the area represented by the boundary node after finite element discretization.
In addition, we apply viscoelastic artificial boundary elements to the truncation boundaries of the solid domain to simulate the wave radiation of a semi-infinite seabed, as shown in Figure 3. e viscoelastic artificial boundary elements are essentially a layer of elements that extend outside the solid domain of the model, and their material properties   Shock and Vibration are specified to realize the function of wave absorption. Under 3D conditions, the equivalent shear modulus G, Young's modulus E, Poisson's ratio μ, and material damping coefficient η of the viscoelastic artificial boundary elements are given as follows [15]: where h is the thickness of the viscoelastic artificial boundary element; R S is the distance from the wave source to the cut off boundary in the solid domain; G and ρ are the shear modulus and mass density of internal solid material, respectively; c P and c S are the velocities of compression and shear waves, respectively; α N and α T are the artificial boundary coefficients, with recommended values in 3D space of 1.33 and 0.67, respectively; and α � α N /α T . Under 2D condition [22], the computational formulas of equivalent shear modulus E, Young's modulus E, and Poisson's ratio μ are all identical to those given in equation (2). Only the material damping coefficient μ holds a different formula, given as follows: Here, the recommended values of artificial boundary coefficients α N and α T in the 2D condition are 1 and 0.5, respectively.

Seismic Wave Input Method.
In the seismic response analysis of an engineering site, commonly, a finite region is intercepted from the actual half-space physical system and the cutoff boundaries are manipulated by using artificial boundaries. en, how to reasonably input seismic waves into the numerical model without affecting the absorption effect of the artificial boundaries becomes an important problem to be considered. Current approaches for solving this problem generally transform incident seismic waves into equivalent seismic loads, thereby inputting seismic energy into the near-field calculation domain. However, traditional wave input methods such as the domain reduction method (DRM) [23,24] and the wave method [25] generally require solving a series of complex formulas to obtain the equivalent seismic loads, which will lead to extreme implementation complexity. Liu et al. recently proposed a new seismic wave input method based on the substructure of artificial boundaries, which is also known as the boundary substructure method (BSM) [16]. Zhang et al. also proposed a substructure method for seismic wave input considering structure-water-sediment-rock interaction [26,27]. e BSM converts the calculation of equivalent seismic loads into the dynamic analysis of a substructure of artificial boundaries, thereby considerably simplifying the calculation processes while ensuring calculation accuracy, especially for seismic wave input in large-scale 3D engineering site models. In this study, we adopt BSM for inputting seismic waves into the model of Zhubi Reef. e calculation processes are   Shock and Vibration shown in Figure 4, and the implementation procedures are given as follows: (1) Establish the reef model and then delete all elements except those connected with the artificial boundary nodes, thereby obtaining a substructure of artificial boundaries. For a given time history of seismic waves, the distribution of the free-wave field can be obtained by conducting a free-wave field analysis according to traveling wave theory. en, by applying the free-wave field motions on the artificial boundary nodes and the adjacent internal nodes, as shown in Figure 4, a dynamic analysis of the substructure model can be conducted. e reaction forces on the artificial boundary nodes can be obtained, which are the equivalent seismic loads.
(2) Apply the equivalent seismic loads on the corresponding positions of artificial boundary nodes in the original reef-seawater model and conduct a dynamic analysis. en, the seismic response of the 3D reef engineering site can be obtained.
A pulse wave shown in Figure 5 and three artificial seismic waves (denoted ad seismic waves A, B, and C) generated based on the bedrock conditions present in the South China Sea, as shown in Figure 6, are selected as the incident waves and are input into the 2D and 3D reef-seawater models through BSM. e acceleration amplitudes of the artificial seismic waves are all set to 1 m/s 2 . Notably, the SV waves input in the 3D model are assumed to vibrate along the short axis of the reef flat (southwest-northeast direction).

Wave Propagation Laws in Zhubi Reef Model
Displacement wave field in the 3D Zhubi Reef model under a vertically incident pulse wave is shown in Figure 7. Meanwhile, we compare the horizontal displacement waveforms on the reef flats calculated by the 2D and 3D models with the normalized parameter x/L as the vertical coordinate, where L is the length of the reef flat. e results are shown in Figure 8. Figure 7 shows that during the early stage of wave propagation, the waveforms in the model of Zhubi Reef are mainly composed of incident SV waves. e first arrival wave reaches the reef flat at approximately 0.6 s, resulting in a significant uniform response over the reef flat. After that, strong seismic responses can be observed in succession around the edge areas of the reef flat. Seismic waves interact with irregular reef topographies and layered reef structures, generating complex scattered waves. ese scattered waves are eventually absorbed by the artificial boundaries, and the reef model tends to be restored stationarily.
e comparison results of the displacement waveforms shown in Figure 8 indicate that the peak displacements calculated by the 3D and 2D models, which are both mainly controlled by the first arrival SV waves, are relatively close. Due to the influence of the local topographies and layer reef structures, however, obvious subsequent waves can be observed at different locations, leading to the extension of the duration of ground motion. Compared with the calculation results of the 2D models, the amplitudes of subsequent waves on Slices 1 and 2 calculated by the 3D models are larger, while the results on Slice 6 are slightly smaller. For the area in the middle region of the reef flat, such as Slices 3 and 4, the displacement waveforms calculated by the 2D and 3D models present a nice consistency.

Amplification of Ground Motions.
To evaluate the seismic magnification effect of Zhubi Reef, we defined the amplification ratios of peak displacement, velocity, and acceleration, R u , R v , and R a , respectively, as follows: where u, v, and a are the displacement, velocity, and acceleration, respectively; the variables x and t denote the xcoordinate and time, respectively; the subscript x represents the horizontal seismic response; and the subscript 0 represents the seismic response on a homogeneous land site under the identical incident wave. Figure 9 shows the spatial distributions of R u , R v , and R a along different slices of the reef flat, with the normalized parameter x/L as the horizontal coordinate. In general, the numerical intervals and the distribution laws of peak seismic responses calculated by the 2D and 3D models agree well with each other but both vary with the observation locations and the input waves. e maximum seismic response over the reef flat tends to appear near the edge area, and the spatial distributions of seismic amplification ratios in this area change drastically. e amplification ratio of the peak acceleration, R a , which is universally considered for the seismic design of reef constructions, can reach up to 3-3.5 in the edge area of most intercepted slices. In the middle region of the reef flat, the amplification ratios of the seismic responses are smaller, and their spatial distributions are relatively smooth. An interesting phenomenon can be observed on the northwest side of Slice 6, where sharp decreases in seismic amplification ratios appear in this area. It can be seen from the local map of Zhubi Reef (see Figure 1) that an inward concave exists in this area. Such irregular topography exacerbates the complexity of the seismic wave field and leads to drastic changes in the ground motion field. In addition, the seismic responses on the reef flat vary with each other under different seismic waves. In general, the distributions of peak seismic response under the incident pulse wave are smoother and smaller than those under the incident artificial seismic waves, which are caused by the Shock and Vibration 5 different amplification effect of the layered reef structure under different waves, according to the previous study [18]. In addition, we compare the spatial average values of the amplification ratio of the peak acceleration R a obtained from the 2D and 3D models (R  Table 2. Step1

Shock and Vibration
On Slices 1, 2, 3, and 4, the peak accelerations calculated by the 3D model are larger than those of the 2D models but smaller than those of the 2D models implemented on Slices 5 and 6. e absolute relative deviations |δ| tend to decrease from the edge region to the middle region of the reef flat. Specifically, for Slices 2, 3, and 4, the relative deviations of the average R a calculated by the 3D and 2D models are less than 8%, while on the areas near the edge of the reef flat, such as Slice 1, the average relative deviation reaches 12.20%, and the maximum value can reach 23.67%. In addition, compared with the seismic responses on the reef flat of Yongshu Reef in the South China Sea obtained in our previous study [18], the peak seismic responses on Zhubi Reef are obviously smaller. is phenomenon can be basically attributed to the different topographies of Yongshu Reef and Zhubi Reef. Compared with Zhubi Reef with a pearshaped reef flat, Yongshu Reef presents a narrower and sharper shape, leading to more significant seismic amplification effect and more drastic spatial distribution of seismic wave field.
Considering the previously analyzed spatial distribution laws of seismic responses on the reef flat, conclusions can be drawn that for the seismic design of constructions on the middle region of a reef site, 2D analysis can be conducted, and the calculated results of seismic response can be amplified by 10% to reflect the 3D topographic effect. However, for the seismic safety evaluation of the edge area of reef flats, where the local 3D topographies change drastically, it is necessary to fully consider the specific 3D topographic characteristics and carry out a 3D seismic response analysis of reef-seawater systems.

Acceleration Response Spectrum.
e reef flats in the South China Sea have been adopted as foundations for engineering constructions such as port terminals, airports, and lighthouses. It is important to provide the local response spectra on the reef flat when conducting the seismic design on those reef constructions. erefore, we further investigate the spatial distribution of acceleration response spectra along the reef flat of Zhubi Reef. e calculated acceleration response spectra with the damping factor of 5% are plotted on 3D images as shown in Figure 10, with the period and the normalized parameter x/L as horizontal coordinates. Furthermore, we calculate the average response spectra on the entire reef flat under artificial seismic waves A, B, and C. e results are given in Figure 11.
It can be seen from Figure 10 that the main components of acceleration spectra calculated by the 2D and 3D models are both distributed within the period interval of 0-1 s. e peak value of spectral acceleration on the reef flat can reach up to 1.2 g and basically appear between 0.1 and 0.3 s, which is considerably amplified and shifted towards the low-period direction compared with those of the seismic waves input from the subsea bedrock. In addition, from the middle region to the edge of the reef flat, the peak value of the response spectrum increases and the period corresponding to the peak value gradually increases. On the northwest side of Slice 6, a significant reduction in the peak response spectrum can be observed, which can be attributed to the irregular topographies in this area shown in Figure 1. From the perspective of the whole reef flat of Zhubi Reef, as shown in Figure 11, the average response spectra calculated by the 2D and 3D models present good consistency with each other. However, some differences can be found in the local distributions of the response spectra. Specifically, on Slices 1, 2, 3, 4, and 6, the response spectra obtained through 3D analysis tend to move slightly towards the low-period direction compared with the 2D results. On Slice 5, an inverse relationship can be observed.

Summary and Conclusions
In this study, a typical reef island in the South China Sea, namely, Zhubi Reef, is selected as the research topic. A 3D reef model considering the fluid-solid interaction as well as the wave radiation effect of semi-infinite seabeds and infinite seawater layers is established, and seismic response analysis of Zhubi Reef is conducted. e distributions of the ground motion field and the acceleration response spectra on the reef flat are studied, and the 3D topographic effect on the site seismic response is investigated through a comparison between 2D and 3D analysis. e following conclusions can be drawn according to the results: (1) Under vertically incident seismic waves, the amplification ratio of the peak acceleration can reach up to 3-3.5 in the edge area of the reef flat, and its spatial variation in this area is drastic. In the middle region of the reef flat, the amplification ratios of seismic responses are smaller and their spatial distributions are relatively smooth. (2) e distribution laws of peak seismic responses in the middle region of the reef flat calculated by the 2D and 3D models agree well with each other. However, for the seismic safety evaluation of the edge area of reef flats, where the local 3D topographies change drastically, it is necessary to fully consider the specific 3D topographic characteristics and carry out a 3D seismic response analysis of the reef site. (3) e acceleration response spectra on the reef flat of Zhubi Reef are significantly amplified and shifted towards the low-period direction compared with those of the incident waves from the subsea bedrock. In addition, from the middle region to the edge of the reef flat, the peak value of the response spectrum increases and the period corresponding to the peak value gradually increases.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.