Liquefaction analysis of soil plugs within large diameter monopiles using numerical modelling

Soil plug formation in open-ended piles due to pile driving is a widely studied phenomenon in onshore applications. The narrow diameters of traditional onshore piles ranging from 0.5 to 1.5 m facilitate plug generation during installation and transfer of enhanced confining stresses across the whole body of the soil plug. Offshore wind monopiles with larger diameters and smaller aspect ratios may not enhance confining stress within the soil plug as effectively as their onshore counterparts. Monopiles are currently the most widespread foundation in the offshore wind sector including in seismic areas. Earthquake shaking can lead to accumulation of excess pore pressure and subsequent liquefaction of the soil retained inside the plug. This research investigates the influence of monopile diameter and confining stress on the development of earthquake-induced excess pore pressure within the soil plug using fully coupled time domain Finite Element Analysis. The results presented in this paper reveal lower accumulation of earthquake-induced excess pore pressure in soil plugs subjected to confining stress. However, the influence of additional horizontal stress caused by monopile driving on the generation of excess pore pressure within the soil plug diminishes as the monopile diameter is progressively increased.


r u
Excess pore pressure ratio t Time in seconds U excess Excess pore pressure zz Strain in the Z direction

Introduction
Open-ended driven piles have traditionally constituted an effective means of supporting structures in the onshore environment. Standard pile diameters range between 0.5 and 1.5 m with L/D ratios of 20 and above. The bearing capacity of an open-ended pile is defined by its base resistance and internal and external skin friction. Pile driving mobilises skin friction owing to the downwards movement of the pile relative to the soil. This leads to a redistribution of the soil stress inside the pile known as soil arching (Shelke and Patra 2008). This originates due to the significantly higher compressibility of sand compared to that of the foundation and provides additional internal friction due to increased horizontal stress. If the pile base resistance is exceeded by the sum of internal soil weight and skin friction, a soil plug will form, for example Henke and Grabe (2008). During driving, the pile will penetrate in an unplugged manner until the base resistance is less than the shaft friction and weigh of soil inside the pile. Following this, the pile will penetrate in a plugged manner (Fig. 1).

Fig. 1 Plug formation mechanism within a pile
Pile foundation technology has been exported to the offshore wind sector, but typical pile diameters have increased to around 7-11 m leading to L/D ratios between 3 and 6. This has given rise to a new type of foundation named the monopile, a large tubular structure made of steel. Most existent offshore wind turbines are supported by monopile foundations. Their simple design, adaptability to different environments and ease of manufacturing (Kallehave et al. 2015) make them ideal for large-scale installation. Monopiles are transported to an offshore wind farm location by vessel and driven into the seabed by hammering. This installation method leads to the formation of a soil plug within the monopile (Fig. 2).
The majority of current offshore wind developments are located in low seismicity areas such as Northern Europe, hence seismicity and liquefaction have not been a primary concern in design to date. The offshore wind sector is however rapidly expanding to other areas of the world such as North America and Southeast Asia, much of which are seismically active, Natarajan and Madabhushi (2022), Seong et al. (2022). Earthquake motions contain a wide range of frequencies, some of which [typically 1-5 Hz as per Madabhushi and Haigh (2012)] can lead to accumulation of excess pore pressure U excess in saturated sands. This may in turn reduce the effective stress of soil and can even cause full liquefaction.
In the classic pile diameter range of 0.5-1.5 m, the soil plug is well confined by the pile and will bear high effective confining pressures (p � ) due to locked-in stresses generated during pile driving. As the diameter of the pile increases, the central portions of the soil plug may not have the same magnitude of horizontal stresses. As suggested by Madabhushi and Haigh (2012), the magnitude of U excess generated depends on the effective confining pressure (p � ) . Thus, it could be postulated that soil liquefaction is more likely to occur in larger diameter piles rather than smaller ones. This paper addresses the development of excess pore pressure in the soil plug of monopiles founded in liquefiable soil subjected to earthquake shaking. This research on the generation of U excess build-up was carried out within a range of soil plug geometries via the Finite Element (FE) code SWANDYNE (Chan 1988). It is hypothesised that earthquakeinduced U excess build-up is considerably lower within the soil plug than in the soil surrounding the monopile. This is thought to be caused by both a reduction in strain due to the structural stiffness of the monopile relative to that of the soil plug and the confining stress exerted by the monopile on the soil plug during installation. This phenomenon was previously observed by Esfeh and Kaynia (2020), who assessed the seismic response of monopiles founded in liquefiable soil utilising a 3D Finite Difference Analysis. Their findings suggest that the difference in U excess build-up between the inside and outside the monopile is caused by the kinematic confinement provided by the foundation. The work includes a liquefaction analysis carried out within a 30 m long and 6 m diameter monopile subjected to static load. Whilst no U excess development is observed within the upper half of the monopile during the seismic motion, a progressive increase in U excess with depth takes place in the bottom half. In fact, the results of the analysis reveal full liquefaction at the bottom of the monopile shortly after the onset of the earthquake. However, the effect of horizontal stresses generated during monopile driving was not investigated by Esfeh and Kaynia (2020).
A reduced rate of U excess accumulation during earthquake shaking may prevent full liquefaction occurring within the soil plug. This may have major implications for the Offshore Wind industry, since the soil within the monopile would partially preserve its strength and stiffness and may therefore provide additional support to the OWT when the soil surrounding the monopile liquefies, particularly in terms of vertical settlement.

Finite element modelling
This paper aims to assess the influence of monopile diameter and confining stress within the enclosed soil plug on the development of seismic-induced U excess and subsequent liquefaction. Hence, the entire soil plug is modelled as saturated loose sand, which is likely to liquefy when subjected to earthquake shaking. A more realistic soil profile would have included some densification of sand with depth.
The FE analyses include four different monopile diameters, 3.33 m, 6.66 m, 15 m and 35.5 m that have been modelled using the fully coupled effective stress based code SWAN-DYNE (Chan 1988). The soil plug depth in all cases is taken as 17.5 m. Nevertheless, the 6.66 × 17.5 m soil plug has been taken as an example for discussion since this is the closest to real monopile soil plug sizes. Each soil plug has been discretised by means of a 2D plane strain mesh (Fig. 4). This is a simplified approach since 2D finite element analysis assumes no strain in the plane perpendicular to the paper. In reality, confining stress applied to the external surface of the soil plug will spread radially across its body, meaning that strain will occur in the 3 dimensions of the plug. The largest amount of strain is meant to occur along the middle cross section of the soil plug and has been therefore, selected for the analysis (Fig. 3). The effects of radial stress variation are outside the scope of the present work and will be included in future research. It must be noted that a full 3D analysis of this boundary value problem can take several weeks to run (Gaudio et al. 2023), whereas the computational time required to complete each run of the 2D analysis reported here is about 3-6 h.
It is assumed that the surface of the steel is rusty and therefore, the internal skin friction is high enough to consider that the soil particles and the monopile wall are fully bonded.
Hence, no interface elements were added to simulate interaction between the inner wall surface of the monopile and the soil within. This assumption is supported by the work of Quinteros et al. (2017) who identified friction angles close to 30 degrees for rough steel, in the same order as the friction angle used in the analyses presented in this paper.
Modelling of liquefaction within the soil plug includes 2 computational stages. The first is static equilibrium where both the initial geostatic stress of the soil and the confining stress resulting from monopile installation are modelled. Earthquake-induced dynamic cyclic loading is subsequently added in the second stage. A range of depths have been selected to assess the response of the soil plug to earthquake shaking (Fig. 4). The static and dynamic computational stages are described in more detail in Sects. 3 and 4.

Static equilibrium
The initial stresses of the soil plug before pile installation and the confining stress resulting from monopile installation are defined in this stage. Firstly, geo-static conditions are achieved by using the boundary conditions of fixed displacement in the horizontal direction (X) and free in the vertical direction (Y) at the nodes located along the left and right vertical boundaries of the model (Fig. 4). Displacement at the bottom boundary nodes is fixed in both X and Y directions, as seismic accelerations will be applied at these nodes later, in the dynamic analyses. ′ xx and ′ zz are the lateral effective stresses generated at a given location depending on the vertical stress ′ yy . The mean effective stress is then expressed as: The geo-static stress state of the soil will change due to the installation of the monopile. Confining stress due to monopile driving is modelled by releasing the fixity in the X direction at the side boundaries and applying horizontal forces in opposite directions at the left and right boundary nodes of the model. In general, installation of the monopile will increase the coefficient of earth pressure (K) . This is defined as: Two distinct stress states are considered for each soil plug model. Soil pressure at rest, with K = 0.5, meaning that the horizontal stress ( � xx ) applied to the side boundaries is half the vertical effective stress ( � yy ) , and K = 1.5, where ′ xx is 1.5 times ′ yy . The latter represents the confining stress exerted on the soil plug during monopile driving into the seabed. Computation of the static equilibrium stage of the soil plugs was carried out for a range of K values comprising 0.5, 0.7, 1, 1.5, 2, 2.5. However, only the K = 0.5 and K = 1.5 stress states were selected in this paper for brevity. Thus, a squeezing effect on the soil plug within the monopile is achieved as it is driven into the seabed. Increased ′ xx involves an increase in the mean effective stress (p � ) within the plug following Eq. 1. The Elastic Perfectly Plastic Mohr-Coulomb constitutive model was selected for computation of the initial and confining stresses. The parameters utilised in the static analysis are listed in Table 1 below. Initially ′ xx , ′ yy and Pore Water Pressure (U) are linearly distributed along the 6.66 × 17.5 m soil plug at rest, i.e. K = 0.5 (Fig. 5). On the other hand, ′ xx and ′ yy are significantly higher on the sides towards the bottom of the plug for the added confining stress scenario with K = 1.5 (Fig. 6). Here, ′ xx is effectively transferred across the body of the soil plug, specially throughout its lower part. In addition, the U distribution remains linear since the potential amount of U excess resulting from the application of additional ′ xx is rapidly dissipated in a high permeability material such as sand. It must be pointed out that at the bottom corners of the models some stress concentrations are observed in Figs. 6 and 7. This is due to the fixed horizontal boundary conditions applied in these analyses, which are required for the subsequent dynamic runs where seismic accelerations are applied along the base nodes. However, the extent of these stress concentrations is rather small as seen in Figs. 6 and 7 and does not influence the liquefaction behaviour of the bulk of the soil plug.
The stress state within the 6.66 × 17.5 m soil plug reached after application of confining stress may be representative of that expected inside monopiles with real diameters. However, the ′ xx distribution obtained for other plug sizes substantially differs from that seen in Fig. 6.
Further analysis revealed that as the width of the soil plug increases, less confining stress is transferred across the body of the plug (Fig. 7). This suggests that the effect of confining stress exerted by the monopile on the soil plug during installation would gradually reduce as the diameter of the monopile increases.

Dynamic analysis
For all the dynamic analyses the geo-static stresses obtained during static analyses are used as the starting stress state. Earthquake shaking is reproduced via an input motion applied to the side and bottom boundary nodes of the soil plug in the horizontal stress ( � xx ) direction, at which fixity is imposed to transfer the earthquake shear stress to the soil body. Displacement in the vertical direction remains free on the sides but fixed at the bottom boundary (Fig. 4).
Flow of water between the inside and outside of the monopile is avoided by using side and bottom closed boundaries. Zero Pore Water Pressure is imposed at the top boundary, i.e., the ground surface. Time step and aspect ratio of the elements may lead to numerical dispersion in wave propagation problems and, therefore, should be carefully selected (Semblat and Brioist 1999). The time step chosen in the dynamic analysis was 1 ms in line with the work by Haigh et al. (2005). Their research recommends adopting a time step one twentieth of the time taken by a share wave to propagate across the smallest element in the mesh at most. As regards the aspect ratio, Potts and Zdravkovic (1999) suggest that the long side of the element should not be larger than five times the short side for rectangular elements. This ratio was maintained between two and three in all models developed for this series of analyses.
The selected input motion is a sinusoidal wave with a maximum acceleration of ± 0.2 g and a frequency of 1 Hz. This loading rate is sufficiently rapid to generated accumulation of U excess within the soil (Madabhushi and Haigh 2012). The P-Z III constitutive model developed by Pastor et al., 1990 was selected to reproduce the soil response against earthquake shaking. The constitutive model parameters utilised in the dynamic analysis are listed in Table 2 below. These parameters were used in a wide range of dynamic analyses reported by Madabhushi and Zeng (1998), Haigh et al. (2005), Madabhushi et al. (2008), Madabhushi and Madabhushi (2014); for the study of liquefaction in boundary value problems. All the analyses in this research were carried out using the same parameters and no modifications to these were done between the analyses.
Occurrence of liquefaction in the soil plug is assessed using the excess pore pressure ratio r u , which is defined as: where r u = 1 means full liquefaction and r u = 0 indicates no U excess generation. The r u time history has been computed for various depths within the 6.66 × 17.5 m soil plug in the K = 0.5 (Fig. 8) and K = 1.5 (Fig. 9) stress states. The input motion applied at the base is For soil pressure at rest (K = 0.5), the r u path displays prominent suction cycles at 1.46 m depth immediately after the start of the earthquake (Fig. 8). Conversely, these do not occur in the K = 1.5 case until after 4 s from the onset of the earthquake, indicating a slower rate of U excess generation. Furthermore, r u values recorded in the K = 1.5 stress scenario are significantly lower than those observed in the K = 0.5 (Fig. 9). This behavior is detected at all control depths except for the point located at 1.46 m depth. In a manner comparable to the K = 0.5 stress scenario (Fig. 8), r u reaches values close to 1 when p � = 0 kPa . However, a higher number of cycles are required to reach liquefaction in the K = 1.5 case, where most of the increment in r u is recorded between 4 and 11 s. On the contrary, r u increases from 0 to virtually 1 in the first 6 s of earthquake for the K = 0.5 case. This is due to the contribution of the extra confining stress added in the static equilibrium stage of the K = 1.5 scenario, which provides higher initial values of p ′ .
The impact of confining stress on the evolution of p′ within the 6.66 × 17.5 m soil plug becomes more evident in Fig. 10, which depicts the stress path of soil in the q − p� space at 1.09 m and 5.47 m depth for the K = 0.5 and K = 1.5 scenarios. At 1.09 m depth the earthquake motion is translated into large cycles of deviatoric stress q that produce a substantial amount of positive excess pore pressure U excess . The mean effective stress p′ decreases accordingly and therefore, the stress path progresses towards the Critical State Line (CSL). At this point, the behaviour of sand changes from contractile to dilatant leading to the development of negative U excess . This results in a slight increase in p′ at the peak of the cycles. However, the overall trend for p′ is downward, meaning that soil stiffness reduces and so does its capacity to propagate shear stress. Hence, the amplitude of the deviatoric cycles decreases as the earthquake progresses. This behaviour is observed for both K = 0.5 and K = 1.5 stress states at 1.09 m depth, but with the difference that it takes more cycles before the stress path approaches the CSL in the K = 1.5 case. Regarding the trends observed at 5.47 m depth, the stress path of the K = 1.5 stress state never reaches the CSL but moves down approximately parallel to it. Therefore, full liquefaction is not achieved at this depth, but the mean effective confining stress reduces from 55 to 25 kPa. This phenomenon could be termed as partial liquefaction. Conversely, in the K = 0.5 case, the stress path starts dropping parallel to the CSL at the beginning of the earthquake, to gradually deviate towards the CSL.
In addition, liquefaction is assessed at various instants in time along the earthquake. It is observed that r u = 0 at t = 0s , meaning no U excess generation (Figs. 11, 12, 13). Following the onset of the earthquake, the values of r u gradually increase from bottom to top within the soil plug, suggesting that full liquefaction initiates a shallow depth and progresses towards the base of the plug. This results from both, an increment in vertical effective stress with depth ′ v0 , and the gradual increase of U excess accumulation rate towards  After the first 5 s of the earthquake, U excess starts building up causing limited liquefaction in the topmost part of the 6.66 × 17.5 m soil plug for both, K = 0.5 (Fig. 11) and K = 1.5 (Fig. 12) stress states. In contrast with this, more than half of the 35.5 × 17.5 m soil plug (Fig. 13) is fully liquefied at this same instant in time regardless of the addition of confining stress. The corresponding r u contour plot displays a high concentration of r u values close to 0 at the top centre point of the 35.5 × 17.5 m soil plug (Fig. 13). These may not be representative of the soil behaviour and are caused by an antinode, which is a point where accelerations generated at both sides of the soil plug are cancelled.
By the end of the earthquake (t = 12 s) , the fully liquefied area (r u = 1) of the 6.66 × 17.5 m plug at rest (K = 0.5) (Fig. 11) is significantly more extensive than that of the same plug with added confining stress (K = 1.5) (Fig. 12). Nevertheless, the complete body of the plug does not fully liquefy in any of the stress scenarios. In contrast with this, the 35.5 × 17.5 m plug with added confining stress (K = 1.5) has completely liquefied at t = 12 s (Fig. 13). This provides clear evidence that the influence of confining stress on the amount of volume of soil that could potentially liquefy within the monopile decreases as the diameter increases.
Mitigation of U excess build-up within the monopile due to confining stress is further corroborated in Fig. 14. This compares the amount of U excess accumulated within the 6.66 × 17.5 m soil plug at rest (K = 0.5) to that generated in the added confining stress Fig. 14 U excess time history for the 6.66 × 17.5 m soil plug with K = 0.5 and K = 1.5 at 5.10 m and 8.75 m depth scenario (K = 1.5) at two different depths. It is clearly observed that the amount of U excess accumulated within the soil plug decreases with depth. Similarly, lower U excess is developed in the K = 1.5 case at both 5.10 m and 8.75 m depth, meaning that confining stress reduces U excess accumulation inside the monopile. This is consistent with the contour diagrams depicted in Figs. 11 and 12. The peak U excess is reached at approximately the same time in both K = 0.5 and K = 1.5 scenarios. However, the maximum value of U excess takes place slightly later at 8.75 m depth, meaning that a higher number of cycles are required to reach the same amount of U excess at a shallower depth. This behaviour matches with that displayed in Figs. 8 and 9.

Discussion
This paper analyses the effect of monopile diameter and confining stress on the accumulation of seismically induced U excess within soil plugs. In order to assess these factors, four different monopile geometries (3.33 × 17.5 m, 6.66 × 17.5 m, 15 × 17.5 m and 35.5 × 17.5 m) have been modelled utilising the F.E. code SWANDYNE (Chan 1988). Each soil plug has been assessed both at rest and with added confining stress resulting from monopile installation. The static equilibrium analysis reveals that confining stress is effectively transferred across the body of the soil plug in the narrowest geometries (3.33 × 17.5 m, 6.66 × 17.5 m). However, this does not apply to the 15 × 17.5 m and 35.5 × 17.5 m, where the influence of confining stress applied at the monopile soil plug boundary due to pile driving on the overall ′ xx distribution is limited. This contrast in the distribution of increased confining stress has significant impact on the development of seismic-induced U excess . In fact, full liquefaction of the entire soil plug is identified for the widest geometries regardless of their initial static stress state. Narrower plugs show restricted U excess accumulation, with full liquefaction occurring only towards the topmost section of the plug. This behaviour becomes more evident when extra confining stress due to monopile installation is added to the soil within the monopile. In addition, it takes a larger number of earthquake cycles for the soil to reach its peak r u in the K = 1.5 case. This is due to the high stiffness of the monopile foundation, and the confining stress developed within the soil plug resulting from pile driving into the seabed.
It is noted that the analyses presented in this paper consider the soil within the monopile only. However, it is hypothesised that if the whole monopile inserted in a sufficiently wide extent of soil was modelled, the bottom of the monopile would have constituted a permeable boundary. Earthquake-induced dynamic cyclic loading would have led to substantially higher U excess generation outside the monopile than inside. Given the pressure gradient, water may have flown towards the inside of the monopile through its bottom boundary, leading to a potential increase in U excess within the monopile in the post seismic period. This pore fluid migration is obviously not captured by the models utilised in this research.
The results presented in this paper reveal that the influence of soil confinement is progressively reduced as the diameter of the monopile increases. In fact, soil plugs exceeding a certain diameter, accumulate sufficient U excess to reach liquefaction regardless of the amount of applied confining stress. This is due to the fact that soil within large diameter monopiles is less influenced by the additional lateral stress. Therefore, a lesser amount of confining stress is transferred across the body of the soil plug during monopile installation. This is translated into smaller mean effective stress and higher accumulation of U excess under earthquake shaking.
Nonetheless, liquefaction is minimised in soil plugs with diameters below 6.66 m when subjected to confining stress. This may involve higher bearing capacity due to residual skin friction within the monopile. On the contrary, soil plugs with diameters close to 15 m and higher, would probably reach liquefaction under the same conditions. Thus, the support provided by internal friction would be gradually lost. Current monopile diameters range between 7 and 11 m approximately, meaning that the soil plug of the largest monopile designs would probably develop enough U excess to reach full liquefaction under earthquake loading.

Conclusions
The influence of monopile diameter and confining stress generated due to pile installation on the accumulation of seismic-induced excess pore pressure and subsequent liquefaction within soil plugs were investigated in this paper. The fully coupled effective stress based F.E. code SWANDYNE was utilised to model four plugs with dimensions 3.33 × 17.5 m, 6.66 × 17.5 m, 15 × 17.5 m and 35.5 × 17.5 m. Soil pressure at rest or additional confining stress caused by pile driving were the initial stress states considered for the soil enclosed inside the monopile. A sinusoidal earthquake motion was subsequently applied to the soil plugs in order to produce U excess . The results of the static equilibrium phase indicated that the distribution of confining stress within the plug is highly dependent on the diameter of the monopile. Added confining stress was effectively transferred across the narrowest soil plugs, i.e. 3.33 × 17.5 m and 6.66 × 17.5 m, whereas the increase of horizontal stress ′ xx in larger monopiles was restricted to the regions adjacent to the lateral boundaries only. Such disparity in the stress distribution had a noticeable effect on the development of U excess during earthquake shaking. The main conclusions derived from the dynamic analyses suggest that for two soil plugs of the same diameter, U excess generation is substantially lower within the plug subjected to confining pressure. This is due to the increase in mean effective stress resulting from the added confining stress. However, it is observed that as the monopile diameter increases, the influence of confining stress on the accumulation of U excess becomes less apparent. This being the case, soil plugs within sufficiently large diameter monopiles may accumulate enough U excess to reach liquefaction irrespective of whether confining stress is exerted on the plug. These finding may be relevant for the offshore wind industry given the need to progressively increase the diameter of monopile foundations. Soil retained within larger monopiles will consequently be more likely to liquefy, meaning that internal skin friction will gradually be lost which may lead to a reduction in bearing capacity.
Funding The authors declare that no funds, grants or other financial support were received during the preparation of this manuscript.

Conflict of interest The authors have not disclosed any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.