A simplified procedure for the prediction of liquefaction‐induced settlement of offshore wind turbines supported by suction caisson foundation based on effective stress analyses and an ML‐based group method of data handling

In this research, a series of non‐linear dynamic finite element (FE) effective stress analyses were conducted to analyze the influence of the suction caisson geometry, ground motion intensity and contact pressure caused by the offshore wind turbine (OWT) on the settlement pattern and seismic demand of the OWT's structure on saturated dense sand. The baseline model and the FE procedure were validated using a database of well‐documented centrifuge tests. However, particular attention was given to the calibration campaign based on the measured system response quantities, such as the settlement, acceleration and pore‐pressure time histories. The FE results identified the contact pressure as an important state parameter caused by the OWT's mass; the governing ground‐shaking intensity measures that play a significant role in the derivation of an analytical framework for predicting liquefaction‐induced OWT settlements during major events are the shaking intensity rate (SIR), Arias Intensity ( Ia)${I}_a)$ and spectral acceleration at a period equal to 1 s (T = 1 s). The results revealed that approximating expressions derived using the modified least‐squares method (MLSM) reasonably capture the complex phenomenon of liquefaction‐induced settlement, with some exceptions at lower SIR values. Finally, to obtain the approximating expressions, the database was combined with a machine learning (ML)‐based group method of data handling (GMDH) that appropriately describes the interplay of multiple properties of the foundation soil, structure and seismic events while incorporating the effect of the interaction between the suction caisson, foundation soil, excess pore‐pressure generation and cyclic shear stresses.


NOVELTY
• The seismic behavior of caisson-supported offshore wind turbines on saturated sand was investigated.
• An insight into the key earthquake and state parameters governing liquefaction-induced OWT settlement was provided.
• The complex interplay of the multiple properties of the foundation soil, structure and seismic events was identified for a given baseline model.
• An efficient machine learning (ML) aided soil-structure interaction (SSI) analysis for a wide range of scenarios was presented.

INTRODUCTION
Wind energy has been emerging as a promising source of energy in developed countries in recent years and has been accepted as a more sustainable and ecological source of energy compared to fossil fuels.One of the major contributors to the cost of an offshore wind turbine (OWT) structure is the foundation support, which accounts for around 20%−25% of the total cost of the structure.In many cases, the cost can reach one-third of the total expense; thus, creating more affordable and cheaper solutions is imperative in order to combat CO2 emissions effectively.One suitable solution to this problem is the suction caisson foundation, which has long been utilized in the offshore oil industry.
As the deployment of offshore wind energy rises globally, there is a significant demand for sites that are appropriate for the installation of wind farms where energy can be generated efficiently.2][3][4] Figure 1 shows the major offshore wind farm developments together with the global seismic hazard map in terms of the peak ground acceleration (PGA).
Seismological and geotechnical hazards may threaten the operational integrity of OWTs in seismically active regions; these hazards may be related to seismological activity in the area, the relative distance from a fault rupture, the forward fault-rupture directivity, the geological composition of the site, the presence of liquefiable soil in the ground, submarine landslides, tsunamis, etc. 5 A typical offshore wind farm is comprised of many different components, including a substation, inter-array cables, export cables and a grid connection.Each of these sub-systems should be operational following a seismic event, and therefore a holistic design philosophy is often recommended. 3,6ylonakis and Gazetas 7 shed light on the controversy and existing misconceptions related to the long-recognized beneficial role of soil-structure interaction (SSI) in the seismic performance of structures founded on soft soil; therefore, supposedly conservative simplification and increased safety margins created by neglecting the SSI effect have been suggested in various seismic codes, 8 NEHRP-97 1997. 91][12][13][14] A few provisions have been developed, such as API-RP-2A-LRFD, 15 IEC-61400-1, 16 DNVGL-OS-C201, 17 and ISO, 18 but they still do not incorporate comprehensive and robust specific guidelines and methodologies that consider the track record of the performance of OWTs experiencing seismic effects and the consequences of soil liquefaction.
To elucidate the influence of kinematic interaction on the seismic loads imparted to an OWT, Kaynia 2 presents a rigorous numerical model of monopile-soil interaction in layered soil, wherein both rotational and horizontal responses at the seabed level were computed as functions of frequency.Knowing that the horizontal component of an earthquake is the main cause of liquefaction triggering, Kaynia 19 demonstrated that wind turbines are, in particular, vulnerable to vertical earthquake excitation due to their rather high natural frequencies in the vertical direction.
He et al., 20 investigated the vertical kinematic response of a large-diameter and thin-walled shell monopile embedded in a partially saturated seabed that was subjected to transmitted seismic P-waves using a semi-analytical method.
Most of the previous research has been exclusively related to either the inertial interaction of ground founded OWTs or the kinematic seismic response of foundations 2,21,22 ; previous studies have mostly focused on how the eccentric the rotornacelle-assembly (RNA) is to the tower top and how the rotary inertia due to blades can influence the structural response F I G U R E 1 Locations of planned and submitted offshore wind farms proposals together with a seismic hazard map showing the peak ground acceleration (PGA). 3 the wind turbines, or they have focused on the forces induced on the piles connected to rigid pile caps following ground motion.A few studies have considered the soil-structure interaction problems of grounded OWT systems in liquefiable soils, but many issues are still uncertain 23,24 ; Haddad et al. [25][26][27][28][29][30] In this study, a comprehensive series of nonlinear dynamic finite element (FE) effective stress analyses modeling the skirted circular foundations supporting OWTs with varied embedment ratios subjected to 20 outcropping rock motions were carried out.
The purposes of the work presented in the current study are (a) To investigate the seismic behavior of caisson-supported OWTs on saturated dense sand as well as to provide insight into the key earthquake and state parameters governing liquefaction-induced OWT settlement.(b) To interpret the complex interplay of the multiple properties of the foundation soil, structure and seismic events identified for a given baseline model and to establish a simplified procedure for estimating liquefaction-induced OWT settlement based on FE analyses using a benchmark model.(c) To present an efficient machine learning (ML) aided soil-structure interaction (SSI) analysis for a wide range of scenarios to aid in the development of the methodology mentioned above.The ML-based group method of data handling (GMDH) is used; this method provides salient insights with respect to the shear-induced OWT settlement mechanisms while eliminating the underestimation of the settlement given by the MLSM for the cases considered within the intermediate range of the SIR.A limitation of the analyses is that they ignore/or not comprehensively incorporate the effects of depth of liquefiable soils, varying suction caisson configurations, sand's relative density or the level of shear strain induced in the foundation soils, seabed permeability, presence of non-liquefiable crust layer, etc., into adopted analytical models.

Centrifuge test configuration of baseline model
The centrifuge tests were conducted by Yu et al., 31 and Wang et al., 32 at the laboratory of Case Western Reserve University to evaluate the seismic behavior of a wind turbine with suction caisson foundation resting on a liquefiable soil.
Figure 2 shows the centrifuge test configuration.All units presented in this manuscript are in prototype scale unless stated otherwise.
F I G U R E 2 (A) Prototype model configuration of a centrifuge test. 31(B) Input earthquake applied to the model.The scaled wind turbine model was installed in a sand-box, which was a rigid container, with a length, width and height of 53.3 cm, 24.1 cm, and 17.7 cm, respectively.Additionally, the carrying payload capacity of such a device is up to 182 kg; however, the gravitational acceleration can be modified by up to 200 g for static analysis and 100 g for dynamic analysis because of the spinning of the centrifuge.
Due to the above payload capacity limitation, the geometrical scaling factor was set to N = 50, which means that the centrifuge model test was performed under 50 g of acceleration.The water table level remained unchanged at 1.5 m above the surface of the sand layer to model realistic offshore environment conditions, as illustrated in Figure 2.
The suction caisson model had a diameter and skirt length of 4 and 1.75 m (i.e., denoted as Test 1), respectively.The suction caisson weight was considered to be 10.7 tons.For the sake of simplicity, a concentrated lumped mass of 10.6 tons was placed at the tip of the tower to simulate the weight (W) of the rotor, blades, gearbox and nacelle of the wind turbine, as shown in Figure 2A.
The rigid container was filled by uniform deposits of Toyoura sand with  50 = 0.17 mm, this sand was poured from the level of 80 cm to ensure that the value of relative density was equivalent to 68% consistently (Yu et al., 2013). 33Moreover, the saturation process was conducted by using a de-aired water system and professional vacuum device for at least 24 h to make the simulation as realistic as possible.More details concerning Toyoura sand are reported in Table 1.It is necessary to mention that 1-D earthquake motion was assumed to be applied directly to the prototype, as depicted in Figure 2B.

FEM modeling considerations
As stated earlier, a nonlinear finite-element model is developed to simulate the seismic behavior of OWTs with a caisson foundation by employing the OpenSees (Open System for Earthquake Engineering Simulation) platform. 34) Soil domain and mesh generation As far as the soil domain is concerned, the pressure-dependent multi-yield surface constitutive model (PDMY01) is utilized to simulate saturated sand.In order to properly simulate the soil domain, QuadUP elements from the OpenSees library are selected for the discretization of the soil domain, since they provide the possibility of simulating the generation of pore water pressure in saturated soil.A QuadUP element is a four-node plane strain element suitable for 2D simulations.Each node has 3 degrees of freedom: the first two are for translations along the x and y axes, respectively, and the third is for the generation of fluid pressure.Nonetheless, not only the kind of element but also the size of the elements plays a key role in obtaining accurate results.The maximum dimension of the elements located in the soil domain can be limited to be smaller than the value of one-eighth to one-fifth of the shortest wavelength considered in the analysis. 35This wavelength can be calculated as follows: where   and   denote the lowest shear wave velocity and the predominant frequency of the input motion, respectively.Thus, meshes composed of QuadUP elements with a size of 0.25 m × 0.25 m, which provide sufficient refinement, are considered, as shown in Figure 3.The discretized model area also had a radius of at least five times the diameter of the foundation, and it consisted of 1908 elements and 2068 surrounding nodes to model the soil medium.

(b) Constitutive model
A constitutive model within the framework of multi-surface plasticity 36 is used here to study the behavior of frictional cohesionless soils.In principle, this constitutive model is an extension of an original multi-surface plasticity concept, with flow and hardening rules 37,38 incorporated.The model has been implemented in the OpenSees software. 38A number of conical yield surfaces is employed; in these conical yield surfaces, the uppermost surface is the envelope of the peak shear strength.However, in the yield function, f is defined as follows: where s =  −  corresponds to the deviatoric stress tensor.In addition,   is defined as the difference between p and , which represent the effective mean normal stress and residual shear strength values, respectively.The other parameters of the formulation include the kinematic deviatoric tensor  and M, which is a material parameter that is represented as a function of the soil friction angle, .
The multi-surface model that is implemented in OpenSees 38,39 incorporates shear-induced contraction and dilation features through a non-associated flow rule.This framework was previously shown to consistently capture hysteretic responses in laterally loaded large slender piles. 40The flow rule is defined separately for stress states that are positioned within or beyond the phase transformation (PT) surface.The plasticity framework also adopts a non-associated flow rule which is typically restricted to the volumetric component  ′′ of the plastic strains, with shear-induced contraction/dilation effects formulated as follows: 37 where  1 and   are scalar coefficients 38 that model the rate of contraction and pore-pressure buildup and the octahedral shear strain, which accumulates over the dilation phase, respectively.In addition,  1 and  2 represent rates of volume increase, 41 while the effective stress ratio is defined by (∕ η), with η tracking the stress ratio of the PT surface.
Table 2 provides a summary of the calibrated model parameters for the liquefiable soil deposit used in the centrifuge experiment.Figure 4 compares the numerically simulated resistance of Toyoura sand to liquefaction triggering with the cyclic simple shear (CSS) laboratory test measurements.Liquefaction was assumed to be triggered when the absolute value of the double-amplitude cyclic shear strain reached 5%.In order to model the structural members of the wind turbine as well as the caisson bucket, the beam-column elements available in OpenSees 34 are used.Due to the fact that the structural elements are expected to remain in their elastic regime under severe seismic excitation, the elasticBeamColumn elements were used, considering the high level of the elasticity modulus, as can be seen in Figure 5. Once the definition of the soil and the structure domain was completed, the interface between the structure and soil was modeled using the equalDOF command technique. 40,42As shown in Figure 5, the nodes of the bucket surface (i.e., slave nodes), are restrained to follow the same displacements as the soil (i.e., master nodes).

(d) Boundary conditions and input excitation
Real structures that undergo seismic forces are subjected to a condition of infinite boundaries in which the dynamic waves generated dissipate from the epicenter toward the free field.Nevertheless, this condition cannot be strictly simulated numerically.Thus, appropriate boundaries capable of dissipating the spurious reflections generated by wave propagation are needed.Based on the centrifuge test, the boundaries appear to be fixed; however, it is a fact that such fixed-boundary conditions cannot provide accurate numerical predictions since they are not able to prevent the reflection of the outwardly propagating waves. 35,43Nevertheless, fixed boundaries are initially needed to ensure the confinement of the soil domain during the analysis.In this study, the boundary conditions have been implemented using the multi-stage method, as depicted in Figure 6.Initially, the bottom boundaries of the model domain are restricted for translations along the x and y directions, allowing the generation of pore water pressure.Regarding the lateral boundaries, only displacements along the x direction are restrained.Second, the gravity analysis is performed using the considered set of boundaries.Subsequently, the horizontal nodal reaction forces generated along the lateral and bottom boundaries are recorded from the model.Then, prior to carrying out the dynamic analysis, the horizontal restraints are removed and replaced by horizontal dashpots at every single node of both the lateral and the bottom boundary, 44 as shown in Figure 6.These dashpots can be modeled by the ZeroLength elements as well as viscous materials available in OpenSees. 34Zhang et al., 35 proposed the following required damping coefficient of viscous elements: where  and  represent the tributary area of each node and the bedrock mass density, respectively.The parameter   is the bedrock shear wave velocity.Finally, to perform the nonlinear time-history analysis, the input motion can be applied to the lateral    and bottom    boundaries using equivalent forces related to the velocity time-history (  ) that are calculated based on the work by Joyner 45 and Ayala and Aranda 46 as follows: where   is the velocity time history, A is the tributary area of each node,  is the density of the rock medium and   is the shear wave velocity of the rock medium.
(e) Damping model Rayleigh damping was included in the soil-structure model.The viscous damping matrix which accounts for small-strain frequency damping with respect to the stress-strain response of the cohesionless soils, 47 can be developed proportionally to the mass [M] as well as stiffness [K] matrices as follows: where the coefficients  0 and  1 can be acquired using the damping ratios   and   for the th and th modes of the system, respectively.The implementation of the Rayleigh approach for the soil profiles with a constant damping ratio was proposed by Park and Hashash, 48 where the coefficients  0 and  1 can be determined using the periods   and   of the system as follows: The two above-mentioned coefficients,  0 = 0.0007879 and  1 = 0.12442, were obtained to cover a full range of load frequencies and the natural frequency of the soil-foundation system, resulting in a minimum damping ratio of 5% at the first natural frequency.

Numerical prediction results
After analyzing the numerical model, the pore water pressure ratio and acceleration spectrum as well as the induced settlements are measured for comparison purposes.A comparison of the pore water pressure ratios from the centrifuge test with those obtained from the numerical model is shown in Figure 7.One can infer from Figure 7 that generally the numerical model is able to predict the trend of all the recordings received from the centrifuge test.The comparison of the results in Figure 7 demonstrates that the adopted numerical modeling method can reasonably estimate the pore water pressure ratio of the system.The computed excess pore pressure (EPP) values are higher than those calculated from the experimental results that are observed in the free-field during the period of strongest shaking, that is, between 2 and 8 s.It turns out that, probably due to the mechanism of sedimentation-solidification postulated by Nagase and Ishihara 49 and Adalier, 50 the excess pore pressure dissipation computed from the PDMY01 model occurs more rapidly than it does in actual case.It is, however, understood from the literature that this less-understood mechanism is likely due to suspended lower-most soil particles regaining their inter-granular contact after liquefaction. 51Keeping the above points in mind, the soil model in the numerical code is divided into four domains, with a virtual nonlinear permeability being defined to account for this phenomenon, as depicted in Figure 3, with the permeability values of 1e-4 cm/s, 2e-5 cm/s, 1e-4 cm/s, and 2e-4 cm/s associated with domains 1−4, respectively.Under the given level of excitation, the soil compaction mechanism caused a steadily reducing trend; the rate of this trend almost became negligible after the end of strong shaking (t ≈ 25 ), as shown in Figure 8A.Finally, the suction caisson settled about 0.184 m; this was caused by an average rate of 1 cm/s during strong shaking.
The computed and measured spectral accelerations are compared at key locations with respect to the OWT structure in Figure 8.The de-amplification of the spectral accelerations at the foundation level and free-field at periods < 1  was captured reasonably well by the numerical model.The period of the soil deposition in the centrifuge test was estimated to be approximately 0.72 s, and the spectral acceleration was de-amplified by nearly 48%.
The soil-structure-induced spectral acceleration was more pronounced at lower periods, and it was captured well.The enhanced de-amplification pattern of soil-caisson-induced spectral acceleration underneath the caisson lid at lower periods can most likely be attributed to the base isolation effect caused by the inclusion of skirts and larger net excess pore water pressure that made the caisson more resistant to rotation; however, it resulted in a lower spectral acceleration.

STRESS-STRAIN RELATION IN BASELINE MODEL
During the computational simulations, soil stress-strain responses were analyzed at key locations alongside the caisson centerline below the lid and free-field (Figure 9).These locations are expected to represent zones of different response features, such as zones with high vertical and lateral confinement, with relatively mild shear within the skirts and major shear strains in the free-field foundation soil (particularly due to the occurrence of liquefaction caused by the imparted seismic excitation).

Free field column
The deformations within the free-field soil at depths of 2.5 and 4 m are associated with a gradual loss of shear strength and the cycle-by-cycle accumulation of shear strain, resulting in large permanent shear strains ranging from around 0.6% to 7.75% (Figure 9D,E).The larger strains occur at shallower depths, while the development of strains near the surface (the shear strains accumulate in both the positive and negative directions) is followed by the degradation of the mean effective stress; the loss of the confining stress occurs at shallower depths, and the confining stress degrades to near-zero values, triggering the occurrence of liquefaction and low levels of shear strength, approximately at the level of 1 kPa.It turns out that the accumulation of large shear strains continues until a depth of 2.5 m.At greater depths, a substantial decrease in the shear strain accumulation is noticed, where at a depth of 4 m, the shear strain accumulation reaches its lowest level, resulting in greater strength.Overall, the trend toward strain-softening is immediately observed after a few initial loading cycles.However, after strong shaking ceased, a regaining of strength occurred, with the stabilizing of the shear strain accumulation.

Below the wind turbine structure
It is uniquely observed in Figure 10 that the asymmetric phases of the dilative response present near the foundation edges spontaneously affect the foundation soil shearing resistance positively, causing the trapped soil within the skirts to exhibit , where  ′  is ef f ective stress).It can be seen that the mean effective stress at z = 1 m degraded to some extent but was far from zero.From the shear stress-shear strain relation, it can be observed that shear strains (  ) developed quickly after the initial shaking and accumulated to a level of 0.4% without further progression.
At greater depths, it can be seen that, due to greater confinement, stress paths initiate at higher effective stresses.As shown in Figure 10, at depths of 2.0 m and 2.5 m, the degradation of the effective confining stress reaches a value of approximately 20 kPa, avoiding the occurrence of liquefaction, while the accumulation of shear strains ended at the levels far below the free-filed.This beneficial effect of SSI is not altogether surprising.Regarding liquefaction countermeasures, Kimura et al., 52 and Adalier 53 were among those that realized the containment capability of sheet-piles, with the interfered drainage capacity as a very effective countermeasure for controlling embankment deformations.The skirt enclosure therefore results in the containment of the migration of loose foundation soil towards free-sides.As a result of this containment remediation, the soil stress-strain response remained in a relatively small range, in contrast to the predominantly contractive behavior at all key locations underneath the skirt tip.
The level of migration and the effectiveness of the skirt are primarily controlled by the hydraulic gradients generated, state variables and motion characteristics.In addition, a comparison between EPPs and stress paths clearly shows that the observed confinement within the skirts and below the tip may be inferred from the areas characterized by the presence of a static shear stress ratio  =    ′  along the foundation edge.Therefore, the foundation soil zone with  ≥ 0.37 induced by an effective stress above ∼ 55 kPa corresponds well to the lower tendency of the shear strain development (Figure 10).
The numerical simulations performed are consistent with centrifuge experiments performed by Bertalot and Brennan 54 and Vaid et al., 55 for relatively dense sand, where the apparent link between the static shear stress ratio, confining stress and relative density is accounted for.

INFLUENCE OF THE CHARACTERISTICS OF THE CAISSON GEOMETRY
A suite of analyses is performed to evaluate the impact of the specific foundation geometry and contact pressure on the magnitude of liquefaction-induced OWT settlement.Four suction caisson-supported wind turbine models with varying aspect ratios (i.e., d/D, where d and D represent the skirt length and diameter, respectively) were evaluated with the

Settlements
From the numerical models, it can be immediately recognized that any of the three alterations of the structure, that is, increasing the diameter of the caisson, increasing the length of the skirt or increasing the superstructure mass, resulted in some sort of behavior improvement in terms of the reducing trend of settlements (see Figure 11).Although the enlarging of the caisson diameter dictates that a smaller vertical effective stress will be induced by the structure in model 2, the volume of the confined soil increases.Hence, the ability of the soil to resist liquefaction is also increased, thereby reducing the liquefaction-induced settlement.Furthermore, increasing the foundation diameter contributes to an increase in the drainage path, which mitigates the significant accumulation of volumetric strains due to dissipation and limits the influence of the ratcheting effect.
In Model 3, in which the length of the skirt is augmented, a noticeable reduction of the settlement is observed.Haddad et al., 4 attributed this to the remedial measures created by the skirts, which caused them to behave as an in-ground, perimetrical, base isolator stiff structural wall and minimized the deviatoric soil deformations (i.e.,  − and  − ) under the wind turbine; this reduced the total OWT settlement by approximately 50% with the increased skirt length, which was in excess of 40%.Additionally, the deeper footing mobilized more sustained net excess pore-pressures and reduced the ability of seismically-induced excess pore water pressures within the soil plug to dissipate rapidly by increasing the drainage path.However, this phenomenon was less obvious when the skirt length continued to increase.However, a skirted foundation is expected to reduce the amount of SSI-induced wind turbine ratcheting ( − ).Depending on the properties of the liquefiable soil, structure and ground motion, if there is a sufficient length of the skirt under the foundation, further increasing this length has been proven to have a minimal impact on the kinematic interaction. 4The water barrier induced by the bucket may reduce the shear stresses transmitted to the structure, reducing therefore the tendency for ratcheting to occur.This can principally be better understood once one investigates the development of excess pore-water pressure (  ) at the four key locations (see the Figure 12).It is immediately observed that the influence of the foundation geometry and state variables on the overall behavior trend is noticeable at the locations within the skirt and beneath the foundation level (P1 and P3 in Figure 2), respectively.
F I G U R E 1 3 Spectral accelerations of select earthquake ground motions.
The lower generation of pore water pressure near the foundation in both models 2 and 3 was likely the main cause of the smaller vertical displacement of the wind turbine models, as seen in Figure 11.Similarly, it was shown by centrifuge experiments and by 1 g shaking table tests that two controversial mechanisms control the shear-induced caisson settlement once the skirt length becomes longer: the bearing capacity failure and the confinement of a larger volume of soil, which contributes more significantly to the reduction of the seismic demand and possible tilting. 4,56An increasing bucket diameter has the same consequences.
A heavier structure in Model 4 resulted in slightly larger moments, which intensified the SSI-induced cyclic shear stresses, excess pore pressures and deformations under the OWT system compared to the Model 1, which had the same foundation area.However, the influence of the weight of Rotor-Nacelle Assembly, represented by a top lumped mass (  ), on the settlement mechanisms is known to be site-specific and primarily controlled by the interplay of multiple properties of the foundation soil, structure, and seismic events.), where V and   hold for vertical loading and ultimate foundation capacity, respectively, is reached.By a series of force-controlled pushover analyses, the bearing capacity is determined.For the purpose of the paper, however, the analyses at low vertical load level  < 10% was preferred to comply with the baseline numerical model assumptions.This is especially relevant to conduct this study at this range, as it is a representative state for operational OWTs, as discussed by Larsen et al. 57 A suite of 20 recorded acceleration-time histories from shallow crustal earthquakes are selected; the characteristics of the selected ground motions greatly capture the seismic response of the baseline case and hence play a key role in provoking the deviatoric-induced OWT settlement.

KEY MECHANISM OF SEISMICALLY-INDUCED OWT SETTLEMENTS
All the selected ground motions are obtained from PEER NGA 58 from representative earthquakes recorded at different locations worldwide; their characteristics are summarized in Table 4.The spectral accelerations of the selected ground motions and the corresponding averages of records are shown in Figure 13.In addition, the spectral acceleration of the synthetic record, which was utilized in the centrifuge tests performed by Yu et al., 31 and Wang et al., 32 is included.
It may be noted that the procedure used to implement the input motions follows the steps as explained above: the velocity time history and shear wave velocity are used to obtain the equivalent forces applied to the boundaries and the dashpot coefficient, respectively.The shear wave velocity required for each case, V S30 ( m s ), corresponds to the shear wave velocity at a depth of 30 m, where the shear wave velocity resembles a recording at a bedrock medium with values greater than 600 m/s.The vertical displacements formed throughout the numerical modeling at the connection between the tower and the lid during each of the events are subsequently computed.
Although the volumetric induced liquefaction settlement can be modeled with purely empirical approaches to account for consolidation-induced strains, 59 the liquefaction-induced settlement of suction caissons atop saturated sand deposits has been found to be governed primarily by deviatoric-induced mechanisms, where the effects of bearing capacity failure and SSI-induced ratcheting types of ground deformation are dominant and associated with the contact pressure, width and embedment of the foundation, motion characteristics, etc.This is in accordance with findings reported by Yoshimi and Tokimatsu, 60 Liu and Dobry, 61 Bertalot and Brennan, 54 Dashti et al., 62 and Haddad et al. 4 There are very few analytical procedures in the literature for estimating shear-induced building settlements.For example, Macedo and Bray 63 developed an analytical procedure based on nonlinear dynamic SSI effective stress analyses for estimating liquefaction-induced building settlement.
In this section, the key features of shear-induced SSI deformation mechanisms identified from effective stress analyses of caisson-supported wind turbines founded on uniform sand deposits subjected to varying contact pressures and a wide range of OWT contact pressure and ground motion characteristics are interpreted, and a simplified procedure is presented in the next section.
For the sake of isolating the individual contribution of the contact pressure, the settlement time histories computed for the varying contact pressures after the implementation of the events R4-Kobe and R10-San Fernando are plotted together with their corresponding Arias Intensity-time histories (see Figure 14A,B).
The Arias Intensity   represents the rate of earthquake energy build-up, and it depends on the amplitude, the ground motion and the duration of the ground motion; it can be obtained from the following equation: where g is the acceleration due to gravity,  is the acceleration and T is the time period.It was, however, found that the time history of the settlements follows the shape of the Arias Intensity-time history of the induced motion, which is consistent with findings reported by Dashti et al., 62 for building footings.
The heavier turbines may intensify larger amount of excess pore pressure generations.During the less intense ground motion (R10-San Fernando), the greater the pressure applied to the soil by the structure, the larger the values of the settlements obtained, as shown in Figure 14B.In terms of the sand's cyclic resistance during the San Fernando event, the resistance to pore pressure generation tends to increase under a higher confining pressure for the   value considered.As a result, the value of the vertical displacements seems to reach a threshold after a certain value of the contact pressure is reached since the energy of motion is not sufficient to counteract the remedial effect of the confining pressure.In contrast to the displacement pattern observed during the moderate San Fernando event, the Kobe shaking was strong enough to compensate for the confining pressure (Figure 14A), thereby intensifying the rocking effect and SSI-induced deviatoric deformations.In this scenario, rather than using the conventionally used parameter, the cyclic stress ratio (CSR), to describe the seismic demand, the authors observe that the OWT settlement is directly proportional to the shaking intensity rate (SIR).This is ascribed to the clear correspondence between the ground motion intensity and its duration and frequency content and the SIR term.
The rate of the Arias Intensity represents the rate of energy build-up.This rate can be quantified by the SIR, which is calculated as follows:  =  5−75 ∕ 5−75 (10)   where  5−75 denotes the change in the Arias intensity from 5% to 75% of its total value, and  5−75 is the corresponding time duration.In addition, the calculated SIR values for all recorded ground motions can be found in Table 4.

Relation between settlements and intensity measures of shaking
In order to find the possible existence of a threshold of the value of the contact pressure above which the intensifying trend of OWT settlement stabilizes, additional masses and contact pressures including m = 60 tons, 70 tons, 80 tons, 90 tons and 100 tons, which provide contact pressures of Q = 46.83kPa, 54.63 kPa, 62.43 kPa, 70.24 kPa and 78.04 kPa, were utilized for all the input motions considered (Figures 15 and 16), keeping all else unchanged.

DEVELOPMENT OF ANALYTICAL MODELS
In this section, an analytical framework for the calculation of OWT deviatoric settlement is proposed using nonlinear dynamic analysis results.Two numerical methods are employed: the modified least-squares method (LSM), which is used as a regression analysis method, and a GMDH.The former was proposed by Farahani and Akhaveissy 64 while the latter was developed by Ivakhnenko 65 based on an artificial neural network.These analytical frameworks provide useful insights into the complex phenomena governing particularly the deviatoric component of the settlement of the foundation soil, where the offshore foundation is present.Finally, the estimated settlement results are compared with observed settlement results from nonlinear time-history analysis (NTHA) to evaluate the accuracy of the developed functions.

F I G U R E 1 5
Relation between the total value of settlements and the value of the Arias Intensity for each the events and contact pressures considered.

F I G U R E 1 6
Relation between the total value of settlements and the value of Shake Intensity Rate (SIR) for each of the events and contact pressures considered.

First method: Regression analyses and functional forms based of the modified least-squares method
For the cases examined in this study, various trial were considered, and the algebraic expression considered is: where  OWT represents the liquefaction-induced OWT settlement (mm), IM1 is the intensity measure and , , ,  and  are coefficients.Regression analyses were carried out using the algebraic equation above (Equation ( 11)) and considering the PGA as the intensity measure (IM) parameter.Henceforth, the ground motion intensity parameters corresponding to the outcropping rocks at the site are incorporated in the proposed simplified algebraic relationship for predicting the shear-induced liquefaction OWT settlement.The functional form constants were found to be  = 0.0282,  = −0.24,= 0.003.= 0.375, = 0.006 and  ≈ 0 when only linear relationship with the contact pressure is considered.
A comparison of the results obtained from the proposed expression in Equation (11) with those from the nonlinear dynamic analysis is shown in Figure 17.

TA B L E 5
The coefficients of the adopted relationship as in Equation (12).Although the proposed model, which uses the PGA as a ground motion intensity measure, results in relatively good prediction of the liquefaction-induced settlement, as expected, the PGA itself exhibits a somewhat noticeable arbitrary uncertainty, particularly when the PGA > 0.8; thus, it is preferable to use the SIR as an intensity measure parameter instead of the PGA.The SIR incorporate the effects of the interaction between the suction caisson, foundation soil, excess pore-pressure generation and cyclic shear stresses.

Coefficient
Furthermore, the regression analyses indicate that the PGA term can be eliminated once the parameter SIR is included.The final form of the recommended equation for estimating shear-induced OWT settlement is: where the   represent the constant coefficients, which are shown in Table 5.
The comparison of the shear-induced settlement estimated using Equation ( 12) and the OWT settlements computed from the effective stress analyses (Figures 15 and 16) during strong shaking is shown in Figure 18.The model predictions are relatively reasonable according to the comparison between the model predictions and the numerical results, with an R-squared value equal to 0.8483, except for some systematic scatter existing where the liquefaction-induced OWT settlements are estimated for low SIR values (i.e.,  < 0.08).Overall, while the model captures the key trends of OWT settlements during large events, the model yet underestimates the calculated settlements for the cases considered within the intermediate range.It turns out that replacing the IM with the SIR instead of the PGA results in a correct prediction of the reduction of the liquefaction-induced OWT settlements when a certain SIR value is exceeded.However, as the available data points with the SIR range of 1-1.75 are sparse to justify the negative inclination branch of the proposed Equation (12), additional care should be taken when using the developed equation for practical engineering aspects within this range.

Second method: A GMDH type neural network-based model for estimating OWT settlement
The characteristics of the selected ground motions greatly capture the seismic response of the baseline case and hence play a key role in provoking deviatoric-induced OWT settlement.In order to identify the important ground motion parameters largely controlling the settlement pattern, the results of the analyses were further examined; it is, however, acknowledged that the relative importance of these parameters can be better explored through Polynomial model or Volterra series,  which is the basis of the GMDH type of neural network proposed by Ivakhnenko. 65A digital appendix provided with this paper 66 has details of structure of GMDH.Favorable IM parameters include the Shaking Intensity Rate (), Arias Intensity (  ) and the spectral acceleration at a period equal to 1 s (  ( = 1 )).Palacios et al., 67 Dashti et al., 62 and Bray and Macedo 68 emphasized the influence of the values of   and  for a given ground motion on the rate of soil structure break-down, excess pore pressure generation, the seismic demand of the structure and SSI effects.Antoniou et al., 69 have conducted parametric study which correlated the seismic foundation settlement with some dimensionless terms.Considering earthquake characteristics, they showed that the influence of PGA and mean frequency of the same, arbitrary shaped ground motion on the induced settlement is more critical.
In addition, the contact pressure () and the mass of the tower structure () were selected as the state parameters.Furthermore, a summary of the statistical parameters of the inputs is reported in Table 6.
In the current study, after the GMDH model is trained, the final optimal function of liquefaction-induced settlement for OWTs is developed.As shown in Figure 19, the best configuration of the GMDH network, which can robustly predict liquefaction-induced settlement, was determined through many trials and errors.A second-order polynomial is considered for each neuron in all layers.Thus, sixteen sub-functions that need to be determined are developed in all layers, as illustrated in Figure 19.
In the first layer of the GMDH network, there are six quadratic polynomials ( 1 ,  2 ,  3 ,  4 ,  5 ,  6 ); each of these polynomials is obtained based on the combination of two input variables using Equation (A4).

𝑀𝑆𝐸 =
∑  (=0) (  () −   ()) 2  (29) A comparison of the results estimated from the proposed Equation (28) (i.e., the outputs of the GMDH neural network) with the results from the nonlinear dynamic analysis (i.e., the targets of the GMDH neural network) is illustrated in Figure 20.The results indicate that the trained GMDH neural network can predict the results obtained from NTHA well.As can be inferred from Figure 20, the generated GMDH perfectly follows the behavior of all output data.Figure 21 illustrates the obtained error, that is, the difference between each sample and its original value.Furthermore, the values for the MSE and RMSE parameters (i.e., Equations ( 29) and ( 32)) for each set of data can be found in Figure 21.The MSE and RMSE values for all data were 9.6818 and 3.1116, respectively.Histograms of the mean error and the standard deviation error for the training data, test data and all data are shown in Figure 21.
Figure 22 shows a comparison of the obtained settlements and those obtained with Equation (28).In order to evaluate the accuracy of the proposed expression, the abscissa represents the settlement obtained by NTHA, and the vertical axis represents the settlement predicted by Equation (28).The error estimation line (i.e., the 1:1 line, or the ratio of the results calculated using the proposed Equation (28) to the results obtained by NTHA when these results are equal) is also depicted in Figure 22.It is clear from Figure 22 that the results estimated by the GMDH model are very close to the 1:1 line, which indicates that the proposed functional form for liquefaction-induced settlement is reasonably able to predict the target settlement values.The values of the regression coefficient for the training data, the test data and all data are 0.8798, 0.9502 and 0.8904, respectively.
As noted previously, a limitation of these analyses is that they mainly focus on the effect of ground motion characteristics and OWT mass, in which the foundation is resting on dense sand, where they ignore other essential factors, those with either deleterious or beneficial impact, which can be used all together to develop a generalized framework for an inherently complex phenomenon.

SUMMARY AND CONCLUSION
Due to extensive coastline and to comply with global climate change agreements, the offshore wind power sector has been expanding worldwide, particularly in natural disaster-prone regions where seismological hazards may threaten the operational integrity of offshore wind assets.The presence of liquefiable soil in the seabed is of primary concern.
Based on the existing literature, the dominant mechanisms of the settlement of buildings supported by shallow foundations on relatively thin deposits of liquefiable, clean sand can conceptually be classified as deviatoric: they involve (a) bearing capacity failure (ε q -BC )and the soil-structure interaction ratcheting effect (ε qratcheting ).
Approximating expressions for estimating liquefaction-induced caisson-supported OWT settlement are proposed.The development of the analytical framework is based on nonlinear dynamic FE effective stress analyses using the OpenSees FE software program with the model, and an efficient ML-and MLSM-aided SSI analysis.Considering 20 earthquake ground motions, the FE analyses are performed by systematically varying OWT and motion properties while all else remains unchanged such as a baseline model (e.g., subsurface conditions).OWT settlement is shown to be greatly influenced by the RNA mass, which can induce varying contact pressures; however, incorporating this effect into the design in practice is not straightforward and requires a thorough understanding of the complex interplay of structural and seismic event properties.
The comparison of the nonlinear SSI response of the baseline model during the moderate San Fernando event with the settlement following the intense Kobe event (the ground motions caused by a reverse fault and strike-slip fault) allowed a deeper understanding of the complex caisson-soil-structure interaction mechanisms, soil nonlinearity, varying contact pressures and intensity of ground motions.For the moderate event with a PGA of 0.17 g and (I a ) of 0.2 m/s, the value of the vertical displacements is shown to reach a threshold value, below which the settlement tends to increase with the contact pressure.However, the dense foundation soil in the baseline model exhibits a greater resistance to seismically-induced pore pressure generation and strength loss under higher confining pressures above the limit.
The fact that the Kobe ground motion compensates for the favorable impact of the confining pressure, in light of higher PGA and I a values, may be attributed to the way that the intensity and rate of ground shaking intensify the seismic demand, which generates larger excess pore pressures and provokes a highly nonlinear SSI response, leading to more extensive seabed breakdown.
To address foundation liquefaction countermeasures, additional turbine and foundation concepts were investigated in which three alterations were made to the structure, that is, increasing the caisson diameter, skirt length and RNA mass.All of these alterations minimized the ratcheting-induced soil deformations under the OWT that may be interpreted as benefits of nonlinear SSI.

F I G U R E 3
Domain distribution and mesh generation of the soil model.

F I G U R E 5
Illustration of the structural elements and interface between the soil and bucket foundation.
(c) Structure model

F I G U R E 6
Schematic of the boundary conditions.

F I G U R E 7
Numerically computed and experimentally measured excess pore pressure and acceleration-time histories at the base of the model (A) at the foundation tip located at the centerline,  3 , (B) in the free field  4 .

F
I G U R E 8 (A) Numerically computed and experimentally measured (A) average foundation settlement (B and C) comparison of FE-computed 5%-damped spectral acceleration responses to those measured in the centrifuge experiment at locations A1 and A2 together with that of the base input motion.

F I G U R E 9
Stress paths and shear stress-shear strain relation in free-field at depths of 0.5 m, 1.0 m, 2.0 m, 2.5 m, and 4.0 m. significantly lower EPP ratios (   =   ′

F I G U R E 1 0
Stress paths and shear stress-shear strain relation recorded at key locations along centered column below the bucket lid at depths of 1.0 m, 2.0 m, and 2.5 m.

F I G U R E 1 1
Comparison of the time-histories of settlements from different performed analyses.specificationslisted in Table3.Model 1 is baseline numerical model with an aspect ratio of 0.43.The second model has an aspect ratio of 0.29 with an increased diameter of 6 m, while all else remains unchanged.The length of the skirt in model 3 is 1.42 times larger than that in the baseline model, making the aspect ratio equal to 0.62.The last model possesses the same dimension as the baseline model but has a greater superstructure mass (38.8 tons).

F I G U R E 1 2
Representative isochrones of excess pore water pressure with respect to depth (L = 1.75 m): (A) Model 1, (B) Model 2, (C) Model 3 and (D) Model 4.

F I G U R E 1 4
Settlement time history and Arias Intensity time history of different events (D = 4 and L = 1.75 m): (A) R4-Kobe and (B) R10-San Fernando.

F I G U R E 1 7
Comparison of data obtained from NTHA with the data given by Equation (11) versus PGA for different contact pressures (Q).

F I G U R E 1 8
Comparison of data obtained from NTHA with the data calculated from the functional form proposed Equation (12) versus the SIR and Q.TA B L E 6A summary of the statistical parameters of the considered database.

F I G U R E 2 0
Comparison between target and network outputs: (A) training data, (B) test data, and (C) all input data.

F I G U R E 2 1
Histogram of the error of the GMDH output, and histograms of the standard deviation error.F I G U R E 2 2 Comparison of obtained from NTHA and those predicted from the proposed Equation (28) for (A) training data, (B) test data and (C) all input data. 31

33 Parameter Particle shape C u C c Specific gravity D 50 D 10 Max. Void Ratio Min. Void Ratio
TA B L E 1 Properties of Toyoura sand (Yu et al., 2013).

TA B L E 2
Calibrated set of parameters of PDMY 01 material.

I G U R E 4
Comparison of computed and measured cyclic resistance of soil to liquefaction triggering for Toyoura sand.
Properties of the models.
TA B L E 3 Characteristics of the input motions.
TA B L 4 Scheme of the developed GMDH network model.