A theoretical and practical framework based on plasticity theory for the drained behavior of single and multiple shallow footings

With the shift of the offshore wind energy sector to deeper waters, demand for the development of more complex foundation solutions, particularly suction bucket – supported tripod/tetrapod and jacket foundations, has increased. This paper is divided into two main sections. The first part comprises a comprehensive review of the performance of circular surface and shallow foundations under combined loading (VHM), and how it can principally be understood in a theoretical framework of plasticity theory. Examination of the considered data suggested that the general assumption of over‐estimated non‐association degree with constant failure surface parameters and increasing vertical load may require further investigation. This may be attributed to the complex interplay of multiple properties such as stress level, soil strength profile and foundation geometry. The existing data in the literature were also used to provide practical guidance for a successful implementation of the elasto‐plastic constitutive relationships in offshore foundation design. In second part of the paper, the suitability of the non‐associated plasticity formulation for a baseline multi‐pod system in H‐M load space was investigated using three‐dimensional (3D) finite element (FE) analyses and not verified. Furthermore, the failure envelopes and hardening law for caissons with different embedment ratios differed from those recommended in the literature were established. Parametric studies of multi‐caisson foundations revealed that the failure mechanism of multi‐bucket foundations under horizontal loading depended greatly on the bucket spacing. The horizontal bearing capacities increased with the bucket spacing until they reached a threshold. Meanwhile, analyses of the multi‐bucket foundation under moment loading confirmed the occurrence of a push‐pull failure mechanism.

F I G U R E 1 (A) Distribution of foundation concepts for European offshore wind turbines up to 2014. 1 (B) Distribution of substructures installed in 2019. 2
Figure 1B shows the distribution of substructures installed in 2019, which clearly indicates the increasing prominence of jackets.This increased use reflects the capabilities of jackets as multiple-foundation systems to stand in deeper water and withstand larger environmental loads.Suction bucket technology was originally developed for the oil and gas industry.Tjelta 3 reported that few suction anchors were installed in the years 1980−1990, but at time, became the preferred foundation type.
5][6][7][8][9] Single shallow foundations have been studied extensively, with findings suggesting the soil-foundation plastic yielding models, [10][11][12][13][14][15][16][17][18][19][20][21] rather than treating the three-dimensional (3D) stress-distribution problem with conventional horizontal reaction concepts developed for piles. 22,23he behaviors of pile-supported and jacket-supported wind turbines differ clearly, and the placement of suction-bucket foundation elements underneath tripod or jacket structures is a new concept in the offshore wind energy sector.In the current study, the need to address the whole nonlinear soil-structure interaction response using macro-element models based on strain-hardening plasticity theory for circular skirted foundations is revisited.However, as the skirted foundations supporting offshore wind turbines are subjected to vertical and horizontal loads as well as moments, the combination of these loads (V, H, M) determines the shape of the yield surface.This surface describes foundation behavior under pseudo-elastic and extremely large deformations within a unified analytical framework.
This paper is divided into two main sections.In Part 1, constitutive relationships for eccentrically loaded offshore foundations are reviewed comprehensively.Generalized failure criteria for the prediction of the drained capacity of circular surface and shallow foundations under general vertical, horizontal, and moment (VHM) loading are addressed (Figure 2).A framework is presented in which plastic flow predictions based on associated and non-associated flow rules, and comparisons of the performance of failure surface and plastic potential parameters, are made.
The second section of this paper builds on discussions presented in previous section, where the appropriateness of the previously adopted framework for a typical tripod supported on multiple caisson foundations was explored using finite element (FE) analyses.Although a few failure criteria had been adopted previously for the monopod configuration under combined loading, we found very few reports on verified failure envelopes for multi-caisson foundations in the literature.
This section validated the consistency of the entire non-associated plastic flow and explored the answer to the key question, whether the normality condition is satisfied in M-H load space or not.The key findings of the hardening law based on the results of FE analyses of the vertical loading as a function of plastic vertical penetration were presented.
Furthermore, the study investigated the alteration of failure surface parameters from those for single shallow foundations and the group efficiency factors corresponding to uniaxial capacities.However, while deriving the non-dimensional The load transfer mechanisms of a-suction bucket supported jacket structure (left) and a mono-bucket foundation (right).
groups for uniaxial bearing capacities of tripod and tetrapod foundations, it was realized that additional scaling relations are necessary to take into account the geometric configuration (i.e., characterizing the spacing distance and embedment ratio).Thus, this section of paper incorporates the group efficiency factors required to study generic multi-pod foundations so that the challenge of using the widely used failure criterion for multiple caisson foundations under general loading can explicitly be answered.Nonetheless, FE analyses are still necessary where understanding the uniaxial capacities of single foundation is crucial.

Failure envelopes
The failure criterion used widely to describe the yield stresses of circular foundations is: where f is the yield surface;   and  0 are the uniaxial tensile and compressive capacities of the foundation, respectively;  is eccentricity parameter that describes the rotation of ellipse in the radial planes.ℎ 0 and  0 represent the uppermost size of the yield surface along the vertical loading coordinate determined by pure horizontal and moment capacities normalized by  0 , respectively.Note that  =0 and  =0 represent intersection values, with the  = 0 and  = 0 axes denoted hereafter as pure horizontal and moment capacities ( 0 ,  0 , and   ,   ).Gottardi et al., 24 evaluated a program consisting of 29 loading tests for circular footings on dry, dense 14/25 yellow Leighton-Buzzard sand, as reported by Gottardi and Houlsby. 25The tests were grouped according to type, stress, and strain history; and consisted of vertical loading beyond peak capacity, followed by unloading-reloading cycles, swipe, monotonic radial displacement with constant ratios of   and 2  (where R is the foundation radius and w, u, and  represent displacements attributed to V, H, and M loads), and constant V and elastic cycles.The estimated friction angles for triaxial compression and plane strain corresponding to   = 75% (  is relative density) at the stress levels of interest below 100 kPa were 42.3 • and 47.8 • , respectively, according to the Bolton 26 procedure.
As a result, the variable term in Equation (1)  1 was considered as function of V and  0 which yields: In continuation of Model C developed by Martin 27 for combined loadings of spudcan in soft clay, Cassidy 10 described the elasto-plastic deformation behavior of spudcan footings on sand using a work-hardening plasticity model and its full F I G U R E 3 (A) Sign convention for load-displacement responses. 28(B) The cigar-shaped yield surface. 29 integration into a dynamic structural analysis program.To establish the model, hydrodynamic loading on jack-up platforms was calculated by integrating wave forces along the leg from the seabed to the instantaneous free-water level, following the well-known Morison procedure and a "deterministic random" wave theory.The final model for the definition of the 3D load-displacement relationship for three loads (V, H, M) and their corresponding displacements (, , ) is shown in Figure 3 and is given as follows: By simply replacing the curvature factors  1 and  2 to account for different variations in curvature at low and maximum vertical stress levels, respectively, a somewhat better fit to the experimental data was achieved.It was further found that slight alteration of the curvature factors would be held with the load path HD/M, but did not pursue this due to data inadequacy.
On the other hand, elliptically shaped capacity envelopes with cone for undrained soil were proposed by various researchers, such as Murff, 36 Bransby and Randolph, 37 Taiebat and Carter, 38 Gourvenec and Randolph 11 ; Cassidy et al., 39 Zhang et al., 40 and Hung et al., 41 The failure criterion in this case can be described as follows: where factors , , and  control the shape and size of strength envelopes that can be illustrated in terms of the aspect ratio L/D (representing skirt length and diameter, respectively), ∕  (where   is the undrained shear strength at the ground surface), and vertical load level  =   0 , respectively.

Hardening rules
The successful implementation of the elasto-plastic constitutive relationships in a foundation design context requires appropriate consideration of the hardening law.3][44] The strain hardening type is expected to predominantly control the foundation response during plastic behavior with a constant The influence of   on post-peak response prediction, according to Houlsby and Cassidy. 45wer-apex failure envelope location (known here as the shallow foundation tensile capacity) if isotropic hardening is assumed for the steadily expanding envelopes.Equation ( 5) is a proposed generalized algebraic expression for vertical loading as a function of the plastic settlement component, which describes the behavior preceding and following the peak 45 : where   corresponds to the ratio of post-peak vertical loading to the peak bearing capacity (   = ) upon the transition of   to infinity, and   refers to the plastic settlement at the maximum load.The   value represents the initial plastic stiffness.Thus, for substantial post-peak softening behavior, where   tends to approach great values, the   value starts to decrease, approaching zero (Figure 4).When using the proposed approach to describing the behavior upon transition to large   values with substantial post-peak work softening (i.e.,   ∼ 0), it is convenient to rewrite Equation (5) by inserting   ∼ 0, which is consistent with that proposed by Gottardi et al., 24,46 and Gottardi and Houlsby. 25

𝑉 =
1 + However, one should be cautious when using the formula for large penetrations   → ∞, which dictates  → 0 but yields a truly unrealistic value.
Alternatively, attempts have been made to relate the work hardening process to not only the vertical displacement term, but also individual components of displacements, as given by: where  1 = 0.5 and  2 = 0.2 are constants corresponding to the yield envelopes postulated in 31,47 for circular plates resting on loose carbonate soil.Equation (7) indicates that the hardening law may be described using Equation (6), with replacement of the plastic vertical displacement term with To establish a vertical bearing capacity formula for skirted footings, Barari et al., 44 and Ibsen et al., 48 compiled data from vertical load tests conducted on dense (  ∼ 80%) sand deposits.A linear form as a function of the circular surface foundation capacity and embedment ratio was proposed: where d and D are the foundation skirt length and diameter, respectively.Amini Ahidashti et al., 49 suggested that the increasing of the skirt length enhanced the bearing capacity in liquefied soils, resulting in a new depth factor exhibited by post-liquefied soils, determined by: Ibsen et al., 14 found that the division of vertical displacement into elastic and plastic components made the interpretation of vertical loading tests for skirted foundations more plausible.

Associated and non-associated plasticity and plastic potential functions
The distinction between associated and non-associated plasticity in the stress-strain behavior of footings will be made clear with the definition of the plastic potential function.For the non-associated plastic flow, the plastic potential function g(V,H,M) defines the direction of incremental plastic strain vector.
Cassidy 10 proposed a generalized approach to the determination of the plastic potential function for a non-associated plasticity framework, with the definition of two association factors ( ℎ and   ) to be consistent with possible changes in the radial (HM) plane due to increases in ℎ 0 and  0 , respectively, thereby reducing the translational and rotational deformations in relation to the vertical displacement. ( 3 + 4 ) )2

𝑔 =
By taking  ℎ equal to   , Equation ( 10) is reduced to 45 : where  ′ indicates the apex of the potential surface.
The association factors show an increasing trend with the development of plastic displacement during the course of loading. 10Subsequently, hyperbolic functions [Equations ( 12) and ( 13)] were developed for eccentrically loaded circular surface foundations on dense sand, and the term  ′ was introduced to control the rate of dependency on plastic strains.Figure 5 depicts the  ℎ and   alterations with the relative magnitude of translational to vertical deformation (  ∕  ) and the diameter times rotational deformation to the vertical displacement (  ∕  ), respectively.Non-associative plasticity is obviously more pronounced in the H-V load space, and to a slightly lesser extent in the M-V load space.
The former term in subscript of  refers to H or M loading coordinate, and the latter corresponds to that occurring at large normalized deformation ratios.F I G U R E 6 Variation in yield surface parameters with the vertical load for surface foundations, according to Model C. 30 After performing a set of regression analyses with the dataset given in, 25 Cassidy 10 suggested the parameters  3 = 0.55,  4 = 0.65,  ℎ∞ = 2.5,  ∞ = 2.15, and  ′ = 0.125, which fit well with the measured plastic displacements.
The failure surface parameters ℎ 0 and  0 were found to degrade nonlinearly, approaching the asymptotic values 0.11 and 0.08, respectively, as the V/ 0 peaked (V =  0 ; Figure 6). 30ℎ 0,  and  0,  data, in terms of the embedment ratio, occurring at near zero vertical load (see Figure 6 for a definition), gathered from, 15,33 are collated in Figure 7.The data were insufficient to determine the analytical expressions with greater accuracy for particular soil densities and embedment ratios of interest (mainly because of slight inconsistency in these values), although the peak failure surface parameters seem to keep increasing slightly with the skirt length.
Notably, over-prediction of the degree of non-association will occur if this decaying trend of the failure surface parameters while  ℎ and   retain their hyperbolically increasing trend is not accounted for during modeling.Despite the increased complexity of potential functions in the literature, which involve various factors that alter with the foundation geometry and soil strength profile, Cassidy 10 obtained a mean value of  ℎ =   = 2.05 that fit the measured data within reasonable tolerance.
For a circular surface foundation, the shape of the yield surface predicted by Model C with a set of yield surface and curvature parameters (  1 =  1 = 1,  = −0.2,ℎ 0 = 0.11,  0 = 0.09), along with the potential function with (  < 1) given by Equation ( 11), is illustrated in Figure 8. Having shown the entire non-associated plasticity formulation for a circular foundation in three-dimensions, it is important to denote that few reports in the literature confirm the validity of the normality condition in the radial (HM) plane.
Two potential surfaces with  ′ =  and  ′ >   have similar shapes despite having different intersection points with the V axis, a clear indication that associated flow is the governing plastic flow mechanism (Figure 8).The latter plastic potential coincides with the yield surface at ∕  = 0.5.Hence, the associated flow rule depicted in Figure 8 is derived by the  ′ = 111.5N,by invoking coincidence with respect to the given yield surface.
Ibsen et al., [14][15][16] also supported the existence of associated plastic behavior in the radial plane, regardless of the embedment ratio, with explicit consideration of the accurate calculation of plastic stress-strain relationships for the adopted yield criterion.
The plastic potential surface for surface foundations was also plotted for two different association factor values ( ℎ ≠   , i. e.,  ℎ = 2.5 and   = 2.1) while all else remains unchanged.The two cases shown in Figure 9  are for two different values of  ′ .Thus, according to Equation (10) and Figure 9, the normality condition is no longer satisfied in the radial plane.The formulation provided enhances the complexity of plastic potential surface description, as it involves several parameters that must be calibrated according to the foundation type and soil strength, and special care should be taken when assessing the skirted foundation responses in different loading planes.
As is evident in Figures 8 and 9 and the above equations, the relative importance of  1 ,  2 ,  3 , and  4 f or the yield and potential surfaces depends greatly on the vertical load level that is described by  = ∕ 0 and  = ∕ ′ .However, the  values in the generalized yield criterion influence the size, but not the shape, of the envelopes in the radial plane.

Some remarks on plasticity models for soil-foundation system (SFS) response
Analogous to the critical state in critical soil mechanics, the so-called parallel point has been proposed in strain hardening models of SFS response; its occurrence upon constant vertical loading and deformation causes the footing to rotate or translate laterally (i.e., at a certain transition point between the punching shear and uplifting where the sliding mechanism predominantly occurs).The parallel point, also defined as the peak state of the potential surface in the V-planes, adopted in several experimental studies (e.g., 50 ), where the foundation can rotate or translate laterally but cannot move upward or downward, is also examined herein to track the potential function envelope and especially to determine a suitable mathematical law for the direction of the plastic strain increments at failure.However, the conditions described above can likely be associated only with the failure surface, due to the occurrence of large deformations upon foundation failure. 16,31hrough analytical and experimental modeling, Byrne 31 and Ibsen et al., 15,16 have developed an interpretation methodology of the evaluation of its initiation.As the parallel point illustrates the peak state of the potential surface along the bi-sectional curves, it requires that , ,   , and    = 0. Furthermore, the plastic vertical settlement at this point is presumed to be nearly zero:   = ∕ = 0.
A numerical study performed to determine the  ′ ∕  ratio revealed a linear relationship between the parallel point location and the apex of the potential surface at a given load state, after the normalization of the obtained parallel point values, 10 expressed by: Additionally, empirical estimation of the parallel point location during failure using the following formula is suggested: Equation ( 14) determines the peak of the potential surfaces in the case of the non-associated flow rule.In contrast, assuming the associated flow rule, where  1 =  3 ,  2 =  4 , and   = 1, the resulting  may be alternately quantified as: On preliminary investigation, the expansion or shrinkage of the potential surface apex due to the alteration of association factor   does not appear to pertain notably to the location of the parallel point under a given load state.Equations ( 14) and (15), in turn, showed that the location of the parallel point defined by two empirical terms ( and ) is more in line with the  factors.
In contrast, the data provided in Figure 10 indicate that   decreases as   increases.Two plastic potential envelopes under given load states appear to follow the shape of the yield envelope, with backward shifting of the parallel point but similar intersection points at the origin.

F I G U R E 1 2
The failure criterion, measured failure data, and plastic displacement increments during failure at load eccentricities M/DH = 0.55, 8.7, and 13 with a small vertical load. 52e apparent dependency of the plastic potential peak state on association factors may allow  to be redefined in combination with other parameters as a guide in the development of a more appropriate formulation in the future.
To isolate the contribution of curvature variable  3 to the shape of the potential function,   and  4 were assumed to be in unity; the resulting functions are shown in Figure 11.The  3 value not only deviated the occurrence peak state, but also caused the envelope to intersect the V axis at the origin more sharply when it was >1; the opposite trend is observed in other apex at  ′  .From this comparison, it can be inferred that  tends to increase.Figure 12 presents proof that plastic deformations are controlled predominantly by translational and rotational, rather than vertical, deformations at low vertical loads, where almost all plastic displacement increments for a shallow caisson (d/D = 0.25) in dense sand tend to be oriented toward vertical.Thus, the peak of the plastic potential surface for a less eccentrically loaded skirted footing (i.e., M/DH = 0.55) may be quantified as a near-zero term,   = 0.05, in agreement with that reported by. 31 , however, continues to increase with the embedment ratio, according to. 51n conclusion, although great effort has been made to understand the complex contributions of curvature and non-unity association factors to peak states in plasticity theory for offshore foundations, this matter remains a major drawback of models, and further research is warranted.

F I G U R E 1 3
Plan view of the multi-bucket foundations.

PART 2: A SIMPLIFIED MODEL FOR MULTIPLE SHALLOW FOUNDATION SYSTEMS
For large turbines in intermediate waters (above 30 m), support structures need to be configured using multiple foundation elements to satisfy the Ultimate LS (ULS).To develop a framework for the prediction of drained HM failure envelopes for a group of foundations in polygonal arrangement, the uniaxial drained capacities and the group effect of multiple caisson foundations (i.e., tripods and tetrapods located at vertices of a regular n-sided polygon) were initially investigated by normalizing the bearing capacity of the multi-caisson system relative to the monopod foundation.This investigation is described in this section.
In the next section, a theoretical framework based on HM failure envelopes and hardening rules is presented.These failure envelopes provided that the failure surface and eccentricity parameters ℎ 0 ,  0 , and , respectively, are dependent on the foundation geometry.
Bordon et al. 53 conducted remarkable work on vertical, horizontal, and sway-rocking stiffness.They showed that mutual interaction among polygonal, rigidly connected foundations produced a helping effect, thus influencing its value more than with the simple addition obtained through Equation ( 17): where   is the stiffness of the foundation group, N is number of foundations, and f emphasizes the terms being for isolated foundations.In the current study, the displacement-control method was applied to systematically determine load-deflection response curves for multiple caisson foundations.Results from more than 150 push-over analyses showed that the bearing capacity of multiple foundation platforms is influenced greatly by caisson spacing and embedment depth; it is expected to differ notably from the single foundation solution due to the interaction among underlying individual buckets.A simplified procedure based on the results of these analyses can be used to estimate the group effect factors.
Figure 13 shows the two multi-caisson configurations used for the numerical analysis.The foundations are connected rigidly and modeled as rigid, embedded cylindrical foundations with diameter D and embedment ratios of 0 ≤ ∕ ≤ 1.To investigate the group effects, the diameter of each individual bucket was kept fixed at D = 6.5 m as the baseline model while the embedded length L was varied from 1.625 to 6.5 m and the spacing ratio S/D was varied from 0.75 to 5.00.The spacing length S describes the spacing between the centers of the individual buckets and the center of the structure [denoted henceforth as the load reference point (RP)].The numerical analyses performed for monopod, tripod, and tetrapod foundations are summarized in Table 1.
For the monotonic loading simulations conducted to obtain the horizontal bearing capacity, a prescribed horizontal displacement was applied at the RP (Figure 14); to obtain the moment bearing capacity, a prescribed horizontal displacement was applied 33 m above the ground surface.
The sign conventions and notations used in the bearing capacity analysis, adopted from, 9,53,54 are shown in Figure 14 and Table 2.The bearing capacity factor   for each push-over analysis was calculated to investigate the group effects of the tripod and tetrapod foundations.It was found by normalizing the horizontal bearing capacity  0 by the number of buckets n, the effective sand unit weight  ′ , and the diameter D to the power of 3. Similarly, the moment bearing capacity factor   was calculated by normalizing the moment bearing capacity  0 by the number of buckets n, the effective sand unit weight  ′ , and the diameter D to the power of 4.

Baseline model
Confidence was achieved by performing a class C prediction exercise with the physical modeling results.Class C prediction, originally postulated by Lambe, 55 was implemented in FE software code using the nonlinear constitutive soil model HSsmall and the Python programming language, to determine whether the simulations successfully captured key failure patterns, verified with bearing capacity results recorded at a trusted high-quality centrifuge facility.This approach allows consideration of the complexities of soil and geometric nonlinearities.Hence, the response of a baseline model 7 under a monotonic loading scenario with silty sand of 11 m thickness (in prototype scale) and ≈ 70% relative density (  ) was assessed numerically (Figure 15).The soil was prepared, and water was added to saturate it.Prior to the sensitivity analyses, calibration was performed based on T2 (monotonic loading).F I G U R E 1 6 Finite-element model of the soil-tripod system.

General modeling considerations
3D FE simulations of a tripod foundation system were conducted using PLAXIS 3D software and Python as the programming language.Preliminary analyses were conducted to determine the mesh fineness and model dimensions required to reach sufficient computational efficiency.
Figure 16 shows the FE discretization of the simulated soil-multiple foundation system.The material properties of steel (E = 210 GPa and  = 0.2) were applied for the linear elastic analysis.The HSsmall model 56 was adopted to simulate the behavior of sand.
The Saemangeum sand parameters were calibrated previously using laboratory element tests. 9The resulting  ′  was 36.1 − −39.4 • ; the dilatancy angle was chosen to be consistent with the formulation of Maranha and Neves, 57 ranging from 1 • to 13 • , and the buoyant unit weight was 9.21 kN∕m 3 .The remaining soil parameters were assumed to be identical to those of silty sand and were defined as functions of the relative density (Table 3).The influence of the underlying ML layer was found to be insignificant due to the formation of shallow failure mechanisms in the upper layer; thus, it was not included in further analyses.As installation likely influenced the in-situ soil conditions and caused apparent loosening, a slight reduction in the relative density from 82.7% to ≈75% was considered.
TA B L E 3 HSsmall parameters Saemangeum sand 9 ; the reference stiffness is for reference pressure   = 100 kPa.The FE analyses were performed in several steps, and interactions between the bucket elements and soil were taken into account using  ′ = 22.6 • .

Model domain and mesh discretization
An equivalent diameter (  ) was defined to represent the circumscribed circle of each bucket's outer edge (Figure 17).The vertical model domain factor   was varied while the horizontal model domain factor   was kept constant (  = 3).  began to converge at a value of ∼4 and stress increment ≤0.2 kPa [Figure 18A].  was set at 4 based on the compromise between the negligible influence from the boundaries and low computational costs.Similarly,   was obtained by assuming   = 4. Again, a compromise was found, giving   = 6 [Figure 18B].

Soil-structure interface properties
To validate the critical state friction angle at the soil-bucket interface, the pull-out resistance for a single bucket in tension was calculated and compared with the measured pullout capacity obtained from centrifuge tests.The pull-out capacity was computed as follows 52 : The experimental pull-out test results 7 suggested that the tensile capacity of the suction bucket was 2.53 MN, yielding an interface friction angle of  ′ = 22.6 • , which corresponds to an interface strength reduction factor of 0.48.Caisson detachment from the surrounding soil and slippage along the interface were admissible forms of geometric nonlinearities throughout the analyses.The drained static moment capacity (  = 92.6MNm) of the tripod was estimated based on Figure 19.The bearing capacities were calculated from the load-displacement curves using the tangential intersection method (Figure 19).Two tangential lines were plotted along the initial and latter points of the curve, where the intersection load was determined as the bearing capacity.The rotation angle achieved was   = 0.256 • .

3.1.4
Loading direction for multi-bucket foundations Kim et al., 7 investigated the most critical direction of horizontal loading with eccentricity for tripod foundations, confirming that the loading direction with one bucket to resist the pullout load (the tension bucket scenario) was the most unfavorable for the tripod.They denoted the load direction as  1 , 7 seen in Figure 15 where two buckets on the leeward side resisted the compression load and the bucket on the windward side resisted the pullout load.The critical horizontal loading direction of  1 was applied for the tripod foundation in this study.The most critical loading direction for the tetrapod foundation was also determined using monotonic loading analyses in the numerical analysis.Figure 20 shows the FE results for the tetrapod loaded horizontally with an eccentricity height of 33 m above the RP in two different directions ( = 0 • and  = 45 • ). = 0 • resulted in the lowest moment bearing capacity and was thus considered to be the most critical loading direction.

3.2
Group effects study

Horizontal group efficiency factor
A prescribed horizontal displacement of h = 1.0 m was applied at the RP. Figure 21 illustrates examples of monotonic responses for the tripod and tetrapod with the selected embedment ratio L/D = 1.00.The tripod and tetrapod were studied with different spacing and embedment ratios, while only the embedment ratio was varied for the monopod (Table 1).The horizontal bearing capacities were analyzed using the bearing capacity factor as in Equation ( 19): Figure 25A shows the horizontal bearing capacity factors  () and  () at different S/D and L/D.Both configurations showed an increasing trend of the   with the spacing ratio until a constant value was reached.For low S/D, the horizontal bearing capacity was found relatively similar and not influenced considerably by the embedment ratio L/D.A disparity in   began to develop beyond S/D > 2.0, when variation in embedded skirt length affected the bearing capacity.
This   pattern proves that the horizontal bearing capacity for multi-bucket foundations depends greatly on the kinematic mechanism that develops upon failure (See Figure 22).As   converged to a constant value at a spacing ratio of S/D = 4.0 with both foundation configurations, the horizontal bearing capacity at S/D > 5.0 was not examined.
Figure 22 shows that both foundations coherently moved horizontally and tilted due to the concentration of deformation on the windward side of the lid.The variation in total deformation indicates that the failure mechanism depends on the spacing ratio S/D.
Although a low spacing ratio dictates the triggering of a failure mechanism similar to that triggered in the monopod, the kinematic mechanism consists largely of horizontal translation at a higher ratio (S/D = 5).
Based on the FE analysis results, functional relationships were proposed for the evaluation of horizontal bearing capacity factors for the tripod and tetrapod structures.The Gaussian function obtained using the curve-fitting technique with two terms [Equation ( 20)] described the variation in   based on the spacing ratio S/D most accurately.The fitted constants of the function as functions of given embedment ratio L/D and structure type are summarized in Table 4.
The horizontal group efficiency factor   =  , ()  , () , which represents the ratio of the horizontal bearing capacity factors of the tripod/tetrapod system and the isolated foundation, is shown in Figure 25B.
The efficiency factors indicate that the effectiveness of the multi-caisson foundation system increased with the spacing ratio, in particular at S/D ≥ 3. The embedment ratio L/D did not show significant disparity until the spacing ratio exceeded S/D > 2.0.
Horizontal group efficiency development was also captured by a two-term Gaussian function with two terms, as seen in Equation (20) in which   was replaced with   and the fitted constants listed in Table 4 were used.The appropriateness of the proposed relationship is shown by dashed lines in Figure 25B.

Moment group efficiency factor
The bending moment applied to the bucket foundations was generated by a prescribed horizontal displacement h with a fixed eccentricity height of 33 m above the ground surface, similar to the conditions in the baseline model.The structures were displaced horizontally until a rotation angle of q = 5.5 • was obtained.Figure 23 illustrates examples of monotonic responses for the tripod and tetrapod structures with the selected embedment ratio L/D = 1.00.In a similar manner to that described above, different geometries with different spacing ratios S/D and embedment ratios L/D were examined.The moment bearing capacities were analyzed using the bearing capacity factor in Equation (21).
increased steadily with the spacing ratio S/D (Figure 25C), indicating that the increase in spacing between buckets affects the moment arm of the push-pull mechanisms in favor of the multi-bucket foundations.The group effect patterns under monotonic moment loading differed markedly from those observed under horizontal loading.Figure 24 shows displacement contours for single and multiple caisson foundations at S/D = 1 and 5 (L/D = 1).The deformations and displacement vectors are concentrated around the front and back sides of the lid, reflecting rotational movement.Remarkably, both tripod systems exhibited a rotational push-pull mechanism, in which the overturning moment was transferred by the buckets as a force couple according to the soil displacement field (Figure 24).Under compression, deformation was concentrated around the entire lid surface because the caissons were being pushed into the sand.
This effect is also reflected in the displacement vectors, where the buckets were pulled out of the sand.When a rotational rigid-body deformation was applied to the RP point, the foundations located at the vertices of the regular polygon coherently moved vertically and horizontally and tilted, seemingly counteracting the effect in the case of rocking rigid-body motion, unlike the boosting effect seen with horizontal loading.
For brevity, the helping effect in other foundations on the same side of polygon was cancelled out to include here.These results agree well with those of the FE analyses of the stiffness correction factors' results for similar groups of foundations reported by Bordon et al. 53 A Simple algebraic equation for evaluation of moment bearing capacity factors for the tripod and tetrapod caisson foundations was proposed based on the FE results [Equation (22)].This equation is proposed based on curve-fitting techniques and expressed as a function of the S/D and L/D ratios.
The fitted constants used, depending on the embedded ratio L/D and structure type, are listed in Table 4.
Figure 25D shows the obtained group efficiency factor   =  () () ;  () is calculated by multiplying the efficiency factor   by the normalized capacity of the single foundation  () obtained using the FE model.As expected, the moment capacity efficiency factor   increased nonlinearly with the spacing ratio S/D, although at a decreasing rate, irrespective of the number of foundations with minimum values of 0.75-1.38 and 0.8-1.8obtained for the tripod and tetrapod structures, respectively, when S/D→ 0.75 (0.25 ≤ ∕ ≤ 1).
Uniquely, under a given S/D, the moment capacity was more pronounced as L/D decreased while n increased.This effect has serious design implications, implying that a tetrapod system should be preferred over a tripod system in the case of a larger moment capacity but no increase in the skirt length.In other words, the factor exhibits, roughly speaking, less sensitivity to the embedment ratio, except at S/D≥ 3.
Given a quadratic function for group efficiency factor approximation, the functional form as in Equation ( 22) was proposed in terms of S/D and L/D.The fitted constants obtained via curve fitting to each function are shown in Table 4.
However, the accuracy of the analytical framework is limited to the full combination of the key parameters deemed for.The appropriateness of the approximation formula is shown with dashed lines in Figure 25D.

3.3
Failure envelopes in the HM plane and hardening rules for the baseline model The bearing strength envelopes for a multi-bucket foundation in drained soil were derived using Equation (1) with the insertion of the normalized tensile capacity  0 = − 3 . As previous research has demonstrated that the combined capacity is affected greatly by the lateral earth pressure on the skirt, different failure parameters depending on the embedment ratio and soil strength were needed to describe the foundation capacity.The failure parameters ℎ 0 ,  0 , and  were derived from FE simulations of the baseline foundation with different embedment ratios, and a vertical load of ∕  ≈ 0.5 was applied to each monotonic loading simulation prior to combined loading.The capacity,   , was determined from vertical compression loading analyses with a spacing ratio of S/D = 3.The recommended prescribed vertical displacement,  = 0.1D, 58 was applied to the tripod foundation with different embedment ratios.The obtained vertical compression capacity   and tensile capacity   of the tripod are shown in Table 5.
Using FE modeling, numerous combined monotonic loads for a tripod foundation were assessed to develop the capacity envelope.Probe loading with the constant ℎ∕ was applied to obtain the capacity envelope under HM loads while keeping the vertical load constant.
The tangent-intersection method was again used to determine the horizontal and moment bearing capacities, with the calibrated soil parameters of the HSsmall model (Table 3) and uniaxial capacities described in the preceding section.
From Figure 26, the plastic displacement increments are shown to be perpendicular to the yield surface at failure; thus, associated plastic flow is present in radial plane.Further, Figure 26 depicts the results in H-M space for embedment ratios of 0.25, 0.50, 0.75, and 1.00 with  2 = 1.0.The failure surface shown for each embedment ratio corresponded well with the FE results, revealing that the shape parameter  1 varied with the embedment ratio and that a better fit was obtained by choosing  1 < 1 and calibrating the factor against the numerical results.This approach to the calibration of the parameters ℎ 0 ,  0 , and  was attempted by varying the aspect ratios (Figure 27).

F I G U R E 2 7
Variation of the yield surface parameters (ℎ 0 ,  0 , and ) for a tripod foundation under a small vertical load. values for tripod system in comparison with those reported for strip footings and isolated foundations.
Consistent with the results obtained for strip and mono-bucket foundations in dense sands, 29,33 the ℎ 0 and  0 values were found to increase with the skirt length: 0 = 1.72 Contrarily, the  0 term given in 15,31 at the near zero vertical load level was shown to be slightly dependent on the embedment ratio for dense saturated or dry sand.On the other hand, the ℎ 0 term shows a linearly increasing tendency which is consistent with those given in saturated sand Barari et al., 33 and dry sand except for embedment ratio of 0.66. 31fter the ℎ 0 and  0 were determined from the size of the yield surface at the intersection with the M = 0 and  = 0 axes, the eccentricity parameter was explicitly determined for the envelopes of combined loading.
Figure 27 shows the linear relationship between the embedment ratio and the eccentricity factor , determined from the analyses conducted with silty sand and a small vertical load.This linear trend has a limiting threshold value of −1.
The eccentricity parameter was also compared with the results of FE analyses of skirted foundations on a multi-layered offshore soil deposit 33 and tests conducted with surface footings on loose sand 31 and medium-density sands. 24,28he experiments conducted with surface footings in dense sand yielded a lower bound value of  = −0.06.This value does not diverge significantly from that obtained in the current study with a low embedment ratio and small vertical load but differs significantly from the value of  ∼ −0.2 reported in. 24,28This observation suggests that  is a function of ∕ 0 .
Of note, the decrease in  with increasing embedment (Figure 27) is likely due to the larger extent of settlement, giving a significant boost to the rotation of the yield surface in the radial plane.Remarkably, a smaller increase in the rotation of the yield surface was observed for the multiple foundation systems.
To establish the hardening law according to Equation ( 5), multiple FE push-over analyses were carried out and the   = 2.18, 1.99, 1.88, 1.85, and   = 14,000 corresponding to the embedment ratios 0.25, 0.5, 0.75, and 1 were determined, respectively.Unlike the post-softening behavior reported in the literature for surface footings, the   values were found greater than unity, with considerably higher initial plastic stiffness.This key finding is in excellent agreement with   ∼12,600   measured by Ibsen et al., 14 for a suction bucket with L/D = 0.75.The stiffness terms are acknowledged to vary markedly among different studies due to the dependency of the behavior on the sand dilatancy, foundation geometry (diameter and skirt length), and stress level.
A limitation of these analyses is that they assume multi-pod systems resting on dense silty sand with constant vertical load throughout the FE studies; therefore, additional research is warranted, which can be used all together to develop a generalized framework.

CONCLUDING REMARKS
A thorough understanding of skirted foundations' responses and plasticity formulations is required to engineer efficient and cost-effective offshore wind turbine foundations.In this work, key remarks on work and strain-hardening plasticity formulations for the ULS analysis of circular-surface and shallow foundations under general VHM loading on drained sand are presented.Salient features and the relative importance and interaction of key factors were studied, and practical guidance for the insightful use of existing relationships based on the comparison of the performance of failure surface and plastic potential parameters is provided.
With the increasing complexity of relationships describing plastic potential surfaces due to the need to calibrate several factors, one must be conscious when assessing skirted foundation responses in the HV, HM, and MV loading planes, which depend greatly on the vertical load ratio, embedment ratio, and soil properties.The generalized plastic potential surface represented in Equation (10) dictates the entire non-associated plasticity formulation of a circular foundation in three dimensions; however, this postulate has not been confirmed in current study with respect to the observed response of monopod and multi-pod systems in radial plane.
The assumption of constant yield surface parameters with increasing vertical loading and hyperbolically increasing  ℎ and   resulted in over-estimation of the degree of non-association.The plasticity formulations for shallow foundations in drained sand do not respond in a fully similar manner during H-V and M-V loading, as is well understood from the scatter observed in the association factors recommended in the literature that has been altered with the development of plastic deformations ( ℎ ∕  ∼ 1.16).Model C demonstrates the degree of non-association given by the nonlinear degradation of ℎ 0 and  0 with increasing vertical load, which peaks when ∕ 0 = 0.Although these trends need to be studied further, the data available in 15,33 suggest that they are valid only for surface footings.The slight dependency of the failure surface parameters on the stress level and embedment ratio suggests that the plasticity models could be further improved.
In this paper, the results of the 3D FE analyses were used to evaluate the bearing capacity of multi-pod foundations in silty sand with group efficiency factors for horizontal and moment loadings.The failure envelopes, where the intersection values with  = 0 and  = 0 axes can be found by the bearing capacity factors in Table 4, were established, with the introduction of closed-form solutions to yield surface parameters and hardening laws as functions of the L/D ratio.
The other key outcomes of this study are summarized as follows: • To accurately estimate plastic strain increments at failure, an understanding of the parallel point state, in which the foundation tends to rotate or laterally translate but does not move upward or downward (i.e.,    = 0,    = 0), is useful.Although preliminary studies showed that the location of the parallel point under a given load state was independent of the   and controlled principally by the alteration of the  factors, its shifting trend as a function of   and  3 was explored in the present work.• The horizontal bearing capacity factor for the tripod and tetrapod structures increased with the S/D and L/D ratios, with foundations at the vertices of the polygons able to move horizontally due to rigid connection.The factor was almost identical at S/D < 1.5, regardless of the L/D; this observation was not confirmed for the tetrapod system.The factor reached a maximum value at S/D ≥ 2 (∕ ≤ 0.5); with continuously increasing skirt length, the maximum value tended to occur at S/D ≥ 4.This steadily increasing trend suggests that the behavior of multi-caisson foundation systems is controlled mainly by a counteracting effect in the case of rocking rigid-body motion, in contrast to the boosting effect observed in the case of horizontal loading.• Gaussian and quadratic relationships corresponding to the bearing capacity and group efficiency factors of horizontal and moment loads as functions of the embedment ratio L/D and foundation structure were proposed.Remarkably, the pattern of limiting values for the horizontal loading factors in terms of S/D was shown to diverge when the moment loading was considered, due to the alteration of the kinematic failure mechanism.• The hardening behavior of the tripod system defined by expressing the size of the yield surface as a function of the plastic vertical settlement differed markedly from that of surface foundations.

A C K N O W L E D G M E N T S
The authors gratefully acknowledge the financial support received from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No 818153.This study is partially funded as part of i4Offshore project (Integrated Implementation of Industrial Innovations for Offshore Wind Cost Reduction).

D ATA AVA I L A B I L I T Y S TAT E M E N T
Some or all data, models, or code generated or used during the study are available from the corresponding author by request.

F I G U R E 5
Alteration of the non-association factors  ℎ and   with   ∕  and   ∕  , respectively.

F
I G U R E 7 ℎ 0, and  0,  for the determination of the bearing capacity of skirted circular foundations under combined loading: comparison of laboratory experiment results and Class A1 prediction according to experiments conducted at the Frederikshavn site in a layered soil profile.Dashed lines show the approximation trend.F I G U R E 8 Yield surface and potential function variation defined in Equation (11) with respect to the V' term (  1 =  1 = 1,  = −0.2,ℎ 0 = 0.11,  0 = 0.09,   < 1).

F I G U R E 9
Comparison of the yield surface resulting from Model C and the potential function proposed in Cassidy 10 in the radial plane ( ℎ = 2.5 and   = 2.1).

TA B L E 1 TA B L E 2
Foundation geometries considered in the SSI analyses.Note: A single caisson foundation is defined by S/D = 0. F I G U R E 1 4 3D view of a tripod caisson foundation.D, diameter of bucket; S, spacing distance; L, bucket skirt length; RP, reference point.Summary of notation for loads and displacements.

F I G U R E 1
Illustration of the tripod system setup.D = 6.5 m, L = 8 m, bucket weight = 4.65 MN, distance between bucket centers = 26.85m, lid and skirt thicknesses = 25 mm.7

37 F I G U R E 1 7
Tripod model dimensions relative to the equivalent bucket diameter,   , and skirt length, L.

F I G U E 1 8 F I G U R E 1 9
(A) Model domain factor for the direction of the load (leeward boundary).(B) Model domain factor for the bottom boundary.Comparison of result obtained with the HSsmall model (   = 0.48) and centrifuge experiment T2 measurement.

F I G U R E 2 0F I G U R E 2 1
Moment-rotation response a tetrapod loaded horizontally at 33 m above the RP with load directions  = 0 • , and  = 45 • , (L/D = 1.0 and S/D = 3.0).Load-displacement response curve for the horizontally loaded tripod and tetrapod with L/D = 1.00.

F I G U E 2
Moment-rotation curves for tripod and tetrapod models with embedment ratio L/D = 1.F I G U R E 2 4 Displacement mechanisms under moment loading (L/D = 1).Tripod caisson with (a) S/D = 1 and (b) S/D = 5.

F
I G U R E 2 5 (A) The horizontal bearing capacity factor   for multi-bucket foundations.  at = 0 the bearing capacity factor for the monopod foundation.(B) Horizontal efficiency   of multi-bucket foundations.(C) The moment bearing capacity factor   for multi-bucket foundations.(D) Moment efficiency   of multi-bucket foundations.F I G U R E 2 6 Direction of plastic displacement increments at failure; FE failure values and proposed failure criterion in the radial plane of the first quadrant for a tripod foundation with S/D = 3.00 are also plotted.∕  = 0.5;  = 1;  1 = 0.97, 0.96, 0.92, and 0.84 corresponding to L/D = 0.25, 0.5, 0.75, and 1, respectively;  2 = 1.0.
plotted elastic and plastic displacements   and   in terms of foundation diameter and obtained average unloading-reloading stiffnesses of   = 15 MPa and   = 6 MPa with   ∼ 4000 TA B L E 4 capacity factors   ,   ,   , and   for tripod and tetrapod caisson foundations in dense silty sand.
Vertical compression   tensile capacity   for tripod foundation with D = 6.5 m and S/D = 3.00.