Numerical study on vibratory extraction of offshore wind turbine monopile foundations under sandy seabed condition

The wind industry has experienced a rapid growth in Europe over the last decades. The early generation turbines were designed for a life of 20 – 25 years. Decommissioning of offshore wind turbines is becoming more important since many of these installed assets are approaching their end of lifetime. In this study, using vibratory extraction of monopile foundations instead of current practice of cutting them is investigated numerically. Correct estimation of extraction force helps operators to choose suitable vibro-hammer and vessel (or crane-barge), which leads to reduction of decommissioning costs. A Coupled Eulerian-Lagrangian (CEL) approach of ABAQUS/ Explicit combined with a modified Mohr-Coulomb (MMC) model are used to find the pile shaft resistance during total removal operation under saturated dense sand condition. The MMC model captures the nonlinear pre-peak hardening and post-peak softening of the dense sand which is not modelled by conventional Mohr-Coulomb model. The VUSDFLD subroutine, which is a user-defined framework, has been used to implement the MMC model into CEL analysis. A parametric study is conducted to analyze how the characteristics of the vibro-hammer, such as its frequency, eccentric moment, and the extraction rate influence the results. The present numerical results show that using proper frequency results in reduction of soil resistance to less than 25% of the initial resistance. However, appropriate hammer with enough eccentric moment and suitable extraction rate, are vital to ensure soil degradation. The results show that the proposed methodology is both robust and straight-forward and it has the potential to reduce computational time which is efficient for engineering applications.


Introduction
The end of life of offshore wind turbines (OWTs) is a subject that needs special attention at the beginning of each offshore wind project; otherwise, the economic impacts can go beyond expectations (Topham and McMillan, 2017).The expected lifetime of the OWTs is 20-25 years (Kerkvliet and Polatidis, 2016;Topham et al., 2019).By the end of lifetime, wind farm operators have to choose an option among extension of the asset lifetime (refurbishment), repowering (replacing old turbines with newer and more efficient turbines) or decommissioning (Martinez Luengo and Kolios, 2015).Even by selecting the first two alternatives based on techno-economic feasibility, regulatory procedures and environmental aspects, decommissioning of wind turbines is still inevitable (Topham et al., 2019).
A decommissioning process starts by isolating the turbine from the grid and de-energizing it.Afterwards, the dismantling proceeds by removing blades, nacelle, and tower (Kerkvliet and Polatidis, 2016;Ortegon et al., 2013).The rest of the process depends on whether the wind turbine is floating, or bottom fixed which requires adopting different strategies for the final stage of decommissioning.Only a few offshore wind farms have been decommissioned so far (Salahshour et al., 2022).In the decommissioning of bottom fixed OWTs, removal of the foundation is the challenging part of the operation.Monopile foundation is the dominant type that is used in 81.2 % of OWTs in Europe (Komusanac et al., 2022).The current solution for removal of monopiles is partial removal which is cutting the pile either externally or internally at or below the seabed level.Offshore cutting practice is time demanding since it requires excavation and dredging operations to reach the determined cutting depth.This is done to create a working space for the cutting tool.In addition, there are other issues related to partial removal of monopiles.Depending on the location, the operator of the farm must then observe the site over a period of several years in order to rule out that the pile stubs do not pose a threat to navigation and fishing in the region.The worst-case scenario would be a net getting caught on the head of a stub which can easily capsize the fishing boat.This can bring life threatening accident to the boat crews.Such seabed irregularities also pose an additional challenge for the approval authorities.When planning and tendering for new offshore wind farms and the associated cable routes, an extra attention should be given to these locations since the remaining pile stubs in seabed pose major hazards.New monopiles should be installed at a distance away from the old foundation.Other type of foundations such as gravity-based foundations may get prone to tilting since the steel body in the subsoil would affect the settlement behavior of the soil.Cable routes would also have to bypass the former locations and jack-up ships could only jack up to a limited extent in these areas.
Another method of foundation removal which is more effective regarding the carbon emission reduction and reducing the operation time (and consequently costs) is total removal of monopiles using vibratory hammers.In addition to the conventional application for installation of piles, vibratory hammers can also be used for extraction purpose.This method was successfully used in decommissioning of Lely offshore wind farm (Tsanova, 2016) as shown in Fig. 1.Each monopile of Lely wind farm was 27 m long, 3.7 m in diameter and with a weight of around 84 tons.Using vibratory hammer resulted in a considerable reduction in removal time of the monopiles.It took only 45 min to extract each monopile (Dieseko Group BV, 2017).The capability of the vibratory technique has been proven in different offshore projects for example project OctaKong in Macau (de Neef et al., 2013;Van Dorp et al., 2019), where 8 synchronized vibratory hammers were used to install monopiles 22 m in diameter and 45 m in length, each weighing 600 tons.The authors believe that the vibratory technique will be used as an effective method which is in agreement with sustainable policies for carbon reduction.
Apart from the time considerations associated with using a vibratory hammer in decommissioning operations, the CO 2 reduction aspect holds significant importance.The recovered monopile steel can be recycled as secondary steel, resulting in emissions that are only a small part of those produced during primary steel production.This highlights the potential of this method for future decommissioning projects.
In geotechnical problems, numerical methods have high potential to reduce expense and time compared to in situ tests or experimental study.However, geotechnical problems such as the installation and extraction processes of piles are still difficult to simulate due to significant and complex deformations (Benson, 1992).The subject of pile driving has been studied widely during the past years.Experimental studies have been carried out by different researchers to investigate the behavior of vibratory driven pile including analysis of the friction fatigue (Moriyasu et al., 2018), the role of the frequency and the actual transferred force from vibrator to monopile (Holeyman et al., 2020), drivability and stress distribution in the vicinity of the pile (Stein et al., 2018), and post-installation dynamic behavior of the pile-soil systems related to different pile-driving techniques (Tsetas et al., 2020).Many numerical studies on the simulation of pile driving have been presented in the literature using various numerical methods.Pile driving is an important aspect of this study because the validation of the proposed model is completed using outcomes of a pile driving field measurement test.
Finite Element Method (FEM) has been a strong tool in numerical modeling of the soil-pile interaction.FEM mesh can be defined using different approaches such as Lagrangian, Eulerian or Coupled Eulerian-Lagrangian (CEL).In the Lagrangian formulation, the computational grid movement and deformation follow the material particles (Belytschko et al., 2014).The Lagrangian method requires significant computational power and can be time consuming.Hence, it is more suitable for small deformation problems.A fully Lagrangian description of the soil has been adopted to simulate deformations and accelerations of saturated sand in the vicinity of the vibratory driven pile (Chrisopoulos et al., 2016;Chrisopoulos and Vogelsang, 2019;Staubach and Machaček, 2019).In the Eulerian formulation, the grid is fixed and allows the material to move freely inside the grid (Benson, 1992).Utilizing a Eulerian mesh prevents mesh distortion problems in the models with large deformation.The solution at the end of each step will be transferred to the initial grid.This formulation based on treating path-dependent material behavior and tracking material interfaces brings extra computational efforts that makes analysis of large deformations problems time-demanding.CEL method attempts to capture the advantages of both the Lagrangian and the Eulerian formulation.A comparison between the Lagrangian formulation and the CEL approach has been conducted for pile jacking into fully saturated soil under partially drained conditions (Hamann et al., 2015).Other applications of CEL include driving of open-ended piles into sandy soils and investigation of the effects of soil plugging on the response of piles embedded in single and multiple layers of soils (Ko et al., 2016) and analyzing the pile installation process and its effects on the pile behavior under a subsequent high-cyclic loading (Staubach et al., 2020).
Arbitrary Lagrangian Eulerian (ALE) is another approach to alleviate the mesh distortion issue in classical Lagrangian approach and it has been used to model cone penetration test in sand (Tolooiyan and Gavin, 2011) and to simulate closed-ended displacement-pile installation in a pressurized calibration chamber filled with sand using Mohr-Coulomb material model (Yang et al., 2020).Multi-Material ALE (MMALE) based on hypoplasticity soil model was successfully applied to evaluate the effects of the frequency in the vibratory installation of tubular piles on the neighboring soil (Daryaei et al., 2018).
Material Point Method (MPM) is developed to overcome the shortcomings of the classical FEM where the (Lagrangian) material points carrying the state parameters move through a (Eulerian) background mesh at which the equations of motion are solved (Phuong et al., 2014).MPM is used to simulate the impact hammering installation of a real size monopile, driven in the North Sea, to demonstrate the capability of this method in simulating the installation of monopiles (Galavi et al., 2019).In another work by Galavi et al. (2017), the mechanical behavior of saturated sand exposed to impact and vibratory pile driving using Mohr-Coulomb (MC) was analyzed.However, these simulations were mostly restricted to ideally drained cases or slow driving processes as encountered for jacking of piles.
Numerical study of pile extraction is rarely found in the literature.Heins (2017) conducted numerical simulations to model complete pile removal.They studied a method called 'free strike' which is used to reduce the shaft friction along the pile when an impact is applied to the pile head.Investigation on the extraction of steel castings for large diameter cast-in-place piles from clayey ground was done by Dong et al. (2018).feasibility of vibratory pile extraction.In the onshore phase, an open-ended pipe pile representative of OWT's foundation was installed and then extracted after a setup period of 6 months.In the offshore test phase, a longer pile was installed and then immediately extracted at 3 different locations in the vicinity of each other.
The classical MC soil model has rarely been used for the numerical modeling of pile driving in sandy seabed due to its unrealistic deformations.The MC model is a simple and fast approach with an acceptable level of accuracy that is based on internal friction angle (φ′), dilation angle (ψ), and cohesion of soil.However, the MC model does not account for the characteristics of non-linear stress-strain response of dense sand including pre-peak hardening and post-peak softening that are affected by accumulated plastic shear strain, loading condition, density, and confining pressure (Bolton, 1986;Hsu and Liao, 1998).Therefore, the MC model may result in inaccurate estimation of stresses and forces.Advanced soil models such as hypoplastic constitutive model, NorSand, and SaniSand have been developed to model the behavior of granular soils.However, a significant amount of computational effort along with calibration based on a large number of parameters poses challenges against the application of these approaches in daily engineering practice.In this study, key features from existing non-linear models are adopted through a simplified approach to incorporate the non-linear effects into a modified Mohr-Coulomb (MMC) model (Bolton, 1986;Nematzadeh and Shiri, 2020;Roy et al., 2015).Large deformation analyses are conducted in ABAQUS/Explicit using the CEL approach.The seabed soil model is built as a user-defined subroutine (VUSDFLD) used in ABAQUS to update the shear strength parameters of the sand based on accumulated plastic shear strain, loading condition, density, and confining pressure in every time increment (Nematzadeh and Shiri, 2020).Monopile removal from saturated dense sand is modelled for both cyclic extraction and direct pulling out of the pile.A parametric study is carried out to evaluate the effects of important characteristics of the vibro-hammer such as vibration frequency and eccentric moment.In addition, extraction rate (or velocity), which is controlled by the crane operator is investigated to find out its effects on variations of soil resistance.In this study, the accumulation and dissipation of pore water pressure is not investigated.An undrained condition is simulated to mimic the soil behavior in the sea and the focus has been placed on the calculation of the pile extraction force.The saturated dense sand modelled in this study is exposed to quick dynamic loading.By considering the high frequency of the vibratory pile extraction (more than 20 Hz), there is no sufficient time for the drainage to happen.It was observed that the adopted simplified approach using the MMC model results in high level of accuracy and a lower computational cost compared to the constitutive soil models.

Base of the numerical model
The numerical model of this study is based on the offshore field measurements carried out by van Dorp et al. (2021).The monopile has a conical shape at the top part with a diameter of 4 meters gradually decreases along the 12.185 m of its head.The monopile length is 56.63 m and it has approximate weight of 301 tons.The monopile was driven to a final penetration length of 17 m below seabed.The water depth at the test location was around 24 m and the soil consisted of medium dense to dense sand.CV-640-VLT-U vibro-hammer, which is a combination of two CV-320-VLT-U hammers, was used to install and then extract the monopile.

Finite element method
In numerical analyses using CEL method, the Eulerian material is tracked as it moves through the mesh by computing its Eulerian Volume Fraction (EVF).A percentage is designated to each Eulerian element, which represents the portion of that element filled with a material (Qiu et al., 2011).The Lagrangian elements have the freedom to move through the Eulerian mesh without resistance until they encounter Eulerian elements filled with material (EVF ∕ = 0).Contact between Eulerian and Lagrangian materials is introduced using a general contact that is based on penalty contact method (ABAQUS documentation, 2008).
Fig. 2 shows the numerical model adopted for the CEL simulations.The void elements that are initially empty of material and material-filled elements are illustrated by green and blue, respectively.The height of void region is modelled in a way to be able to capture soil heave during driving and extraction of pile.It is worth mentioning that only elements filled with material can contribute to the global force equilibrium.The soil domain is considered large to ensure that the effect of boundary conditions on soil flow is diminished.
Due to pile negligible deformations compared to soil deformation, the pile is modelled as a rigid body using Lagrangian elements.The bottom boundary of the Eulerian region is constrained in vertical direction while the side wall is constrained in horizontal direction.Kinematic contact formulation is used to model the pile and the subsoil interaction.This formulation suppresses the transmission of the tensile contact stresses and only shear stresses are transmitted, which includes Coulomb's friction law.The friction coefficient is calculated based on μ = tan(δ), where δ is the friction angle between the pile and soil (Helwany, 2007).However, CEL provides better results for cases of smooth (μ = 0) and rough (μ = 1) (Dutta et al., 2015).Thus, a rough pile is considered in this study.To ensure the accuracy of the numerical simulation, initial geostatic stresses balance was introduced in a separate step by applying gravity load to whole soil domain.The focus of the present study is to investigate the pile-soil interaction; hence, the effects of surcharge force due to hydrostatic pressure of sea is neglected in the simulations.

The adopted constitutive soil model
Mohr-Coulomb (MC) material model is more reliable for noncohesive soil since the pore water pressure drains under the loading in this type of soils.While the MC model assumes that plastic strains only occur when the stress state is on the failure (yield) surface, experimental studies demonstrate that they develop much before failure.Once the stress state reaches the yield surface, the soil deforms with a fixed dilation angle, and any changes in stresses inside the yield surface just causes elastic strain.To capture this behavior, several constitutive models have been proposed in the past (Dafalias and Manzari, 2004;Gajo and Muir Wood, 1999;Prevost, 1985).In the present study, non-linear models are adopted to incorporate the effects of non-linear stress-strain behavior that is not considered by the built-in MC model in ABAQUS (Bolton, 1986;Roy et al., 2015;Vermeer and De Borst, 1984).The present model incorporates the important features of stress-strain behavior of dense sand such as the nonlinear pre-peak hardening and post-peak softening of the internal friction and dilation angles with plastic shear strain, loading conditions, density and mean effective stress.A large deformation CEL analysis of ABAQUS/Explicit is unable to model pile driving and extraction using MC material model.However, the adopted model makes performing effective stress analysis possible through implementing the non-linear aspects of a dense sand (Nematzadeh and Shiri, 2020).Moreover, the dense sand in this study has been verified against the published dense sand tests (i.e., test P06) reported by Yang (2009).

Flow rule
Flow rule is a relation that predicts the plastic strain of the material induced by the stress.Yielding in dense sand occurs after the material reaches to the peak friction and dilation angles.Bolton (1986) conducted a series of triaxial tests to investigate the dependence of the peak S. Salahshour et al. friction angle (φ ′ p ) on mean effective stress (p ′ ) and also the relationship between the peak dilation angle (ψ p ), peak friction, and critical friction (φ ′ c ) angles.The following flow rule is used for describing behavior of the dense sand: where I R is the relative density index and I D is the relative density (Dr).Q and R are equal to 10 and 1, respectively, which are the best fits for most of the laboratory tests, but they may vary with the type of sand and p ′ .The subscripts P and C describe the peak and critical states, respectively.
A is a constant equal to 3 for the triaxial condition and 5 for the plane strain condition.φ ′ c is the critical friction angle and is taken as 31 • .K is equal to 0.5 for the triaxial and 0.8 for the plane strain condition (Bishop, 1961).In this study, as the pile is not long enough to be considered as a plane-strain object to capture the complexities and non-uniformities in the actual pile-soil interaction, a three-dimensional version of the flow rule is adopted.

Strain hardening and softening
The following equations are used to approximate nonlinear hardening of the mobilized friction angle as a function of accumulated plastic strain (Vermeer and De Borst, 1984).It is assumed that the mobilized friction angle increases from 0 • to φ ′ p when γ p p is reached.
The parameter φ ′ in is the initial friction angle corresponding to the beginning of the plastic response.The mobilized friction angle during post-peak softening can be expressed as follows (Hsu and Liao, 1998;Vermeer and De Borst, 1984): Fig. 3 illustrates the geometrical form of Equations ( 9)-( 12) for various magnitudes of relative density and effective mean stress of 40 kPa, where the post-peak softening inflection point is located at a shear strain of γ p c / ̅̅̅ 2 √ which is greater than γ p p (Nematzadeh and Shiri, 2020).The initial part of the curves depicts the strain hardening that starts from the pre-yield zone in which φ ′ and ψ increase from φ ′ in and ψ in to φ ′ p and ψ p at γ p p .It is important to note that for dense sands with low confining pressures, the peak strengths are mobilized at relatively low strain levels.Thus, the plastic shear strain that corresponds to the peak friction angle can be expressed as a function of both the relative dilatancy and confining pressure (Hsu and Liao, 1998;Lee and Seed, 1967;Lings and Dietz, 2004;Tatsuoka et al., 1986).By adopting this approach, the following equations are valid to consider the effect of density and stress level on γ p p (Roy et al., 2015): where γ p c is the strain-softening parameter and p ′ a is the atmospheric pressure (100 kPa).The soil parameters m, C 1 and C 2 can be determined by performing a series of triaxial or shear experiments at various confining pressures and densities.
Poisson's ratio (ν) and Young's modulus (E) represent the elastic parameters of the soil.The Poisson's ratio of 0.2 has been considered as a representative value for dense sand (Jefferies and Been, 2015).The Young's modulus is variable through the depth of the soil based on the initial confining pressure and is expressed by the following power function (Janbu, 1963): where k 0 is the at rest lateral coefficient, and K and n are the power function parameters equal to 326 and 0.86, respectively.H is the depth of the soil and γ ′ is the effective unit weight.
A user-defined subroutine (VUSDFLD) is developed and incorporated into explicit CEL analysis to capture the non-linear stress-strain behavior of dense sand.In every time increment, the stress and strain components are fed to the subroutine.The subroutine first calculates p′ using the stress components.Then the principal strain components are obtained using the VSPRINC utility.The gradient of the principal plastic strain is then determined as the difference between its minor and major components.The sum of the Δ γ p during the analysis time yields the γ p .The parameters γ p and p′ are defined as field variables FV1 and FV2, respectively.By using the MMC model equations, the mobilized φ′ and ψ are defined as a function of FV1 and FV2 in a tabular form in the model input file.Afterwards, the ABAQUS software utilizes the information from the subroutine to update the values of φ′ and ψ based on the field variables.This process takes place simultaneously during the analysis and includes the non-linear stress-strain behavior of the dense sand material (Nematzadeh and Shiri, 2020).

Verification of numerical model performance
The practicability of the soil model was investigated for an ice gauging event in dense sand (Nematzadeh and Shiri, 2020).The results of the sand model with two different relative densities (50.8% and 39.0%) were compared with the results of free field ice gouging testing of Pipeline Ice Risk Assessment and Mitigation (PIRAM) Joint Industry Program (JIP) (Yang, 2009).Moreover, three parameters of MMC model for dense sand were calibrated through a set of triaxial or simple shear tests at different confining pressures and densities (Roy et al., 2015).The sand parameters of the MMC are given in Table 1 and will be used to represent dense sand behavior.
The MMC model has been able to capture the cyclic oscillations of the reaction forces after entering of the ice to the seabed which were observed in experimental studies.Additionally, to further investigate the ability of this model for characterizing sand behavior under cyclic loading, a validation with field measurements for driving of monopile is carried out.In the test conducted by van Dorp et al. (2021), the 56.63 m long monopile was driven 17 m into the seabed and then extracted using a vibro-hammer which its characteristics are provided in Table 2.
The used hammer in field measurement tests can produce a maximum centrifugal force of 13.756 MN when it works with the highest frequency.However, during the driving phase, the hammer frequency was changing between 19.5 Hz and 23.3 HZ.
In the present study, the displacement-controlled technique is employed since the primary objective is to determine the amount of force necessary to extract the pile.This methodology is also adopted in validation of the model for pile driving case.The pile displacement, both during driving (for validation purpose) and extraction (for analysis purpose) comprises two components, i.e., the cyclic oscillations and the rate of driving/extraction.Field measurements during the pile driving reported by van Dorp et al. ( 2021) indicated that it took approximately 7 min to install the monopile to a depth of 17 m in the dense seabed sand.It took almost the same time to retrieve the monopile during the extraction phase.Thus, a velocity of 40.5 mm/s is employed for pile movement in the present numerical simulations.In order to apply the cyclic oscillations induced by the vibro-hammer to the pile, theoretical formulas from Holeyman (2002) are used which is elaborated here.The maximum centrifugal force F c created by the vibro-hammer depends on the eccentric moment of two (or another even number) of contra-rotating eccentric masses and the frequency f d and is obtained as follows (Holeyman, 2002): where M e and ω are eccentric moment and angular frequency that are calculated as follows: m e,i is the mass of counter-rotating masses and r e,i is the distance between the center of gravity of individual mass and the axis of rotation.
The theoretical peak displacement amplitude of pile and vibro-hammer can be calculated as follows: The term m dyn refers to the mass of all components that vibrate, such as the pile and the vibro-hammer, and is therefore known as the dynamic mass.This cyclic oscillations with amplitude of s 0 is applied to the pile in addition to the rate to create cyclic driving or extraction.Assuming that the amplitude of the pile's vibration is equal to the value obtained from Equation ( 19) is overly optimistic since soil typically exerts a damping effect (Jonker, 1987).Fig. 4 shows the present numerical results of pile driving as compared to the field measurement data for the purpose of validation.The amount of force needed to drive the pile to the particular depth of 17 m is obtained from numerical results and is provided in Table 3.This force is the sum of the inner and outer shaft resistances, in  addition to the resistance at the tip of the pile.The pile driving simulation does not result in any soil plugging and only small non-uniform heaves are observed inside and outside the pile on the top of the soil.This indicates that the impact of soil displacement caused by the pile's penetration can be disregarded in the boundary condition for the extraction phase.Thus, the state of wished-in-placed (WIP) is considered for extraction simulations.

Results and discussion
This section shows the results of the present numerical study.In the first part, a comparison is made between cyclic extraction and direct pulling out of the monopile i.e., without using vibro-hammer.Then, in the second part, a parametric study is carried out to analyze the effects of important vibro-hammer parameters such as hammer generated frequency, eccentric moment, and pile extraction rate on soil resistance.

Cyclic extraction versus pulling out
In order to extract the pile from the seabed, it is required to overcome pile axial shaft resistance, weight of vibro-hammer, submerged weight of monopile, marine growth, and pile plug if plugging has happened.Even if a plug has not been formed during installation phase, there is still possibility of plug formation during the operational phase of the monopile, because of the absence of inertia of the soil column during this phase (Fattah and Al-Soudani, 2014).This implies that the static loading due to soil set-up (or settlement) can result in the plug formation.In the present numerical simulations for dense sand, no plugging is observed during driving and extraction and the reason can be related to the large diameter of the monopile.Only the soil resistance is calculated in the present simulations.Hence, only this parameter is compared to the analytical formula found in the literature (API, 2007;DNV, 2014;ISO19902, 2007).The weight of the vibro-hammer, the submerged weight of monopile and the marine growth effect are ignored.
The pile shaft resistance under tension, Q f ,t , represents the frictional force that soil exerts on both inner and outer surfaces of the pile (see Fig. 5).It consists of the area of the pile multiplied by the soil's unit skin friction, which is dependent on the soil type, the submerged unit weight of the soil and the penetration depth.There is no distinction between the inside unit skin friction, f i,t , and the outside unit skin friction, f o,t .The difference between the pile inner and outer shaft resistances is just due to the surface area, i.e., depending on D i and D o (Vugts, 2016).This is shown in the following equations, The unit skin friction of pile in sand is determined according to the following equation (API, 2007;Vugts, 2016):

Table 3
Soil resistance at penetration depth of 17 m.

Allnamics prediction
The present study (ABAQUS)
The unit skin friction resistance originates from the horizontal effective stress, K 0 σ ′ v (z), and Coulomb friction, tan(δ), or, in other terms, the horizontal effective stress multiplied by value of β (ISO19902, 2007).At large depths, the skin friction no longer increases, and the value is constrained to f lim (Vugts, 2016).Although the exact explanation for this phenomenon is not fully comprehended, it could be associated with the fact that soil particles, under extremely high levels of stress, may shift and reorganize themselves in an attempt to alleviate the stress.The values of limiting stress and β are taken from several design standards and are summarized in Table 4.
For an unplugged scenario in sandy soil, the total shaft resistance of a pile, which is the sum of the inner and outer shaft resistances, can be computed using Equation ( 22).The depth at which the limit skin friction is reached, z f,lim , and the shaft resistance thus increases linearly, lies between 20.5 and 23.0 meters for almost all sand types.
Based on the theoretical method, the total pile shaft resistance in dense sand for a pile with 4 m in diameter and 17 m embedded length is shown in Fig. 6 for two different cases (This pile is considered as the reference pile in this study).The curve with lower values is drawn based on DNV's recommendations where the shaft outer resistance is multiplied by a coefficient equal to 0.625 (DNV, 2014;Meijer, 2022).
Soil resistance can be reduced by using vibration, and it can be estimated using the β method (Jonker, 1987) analytically as outlined in Equation ( 23).The β factor is the ratio of the residual strength after cyclic loading to the initial strength, as displayed in Fig. 7.The vertical axis of Fig. 7 represents the soil strength, and the horizontal axis represents the displacement.This factor is dependent on the soil type.The coefficients β o and β i are the ratio of soil resistance degradation outside and inside of the pile, respectively.Jonker (1987) provided some indicative values for different types of soil for pile driving.The value of β may be differ for extraction condition.Based on the post analysis of vibratory driving and extraction records, a value between 0.2 and 0.25 is suitable for fine to coarse-grain sand.Consequently, for the considered reference monopile, the value of pile shaft resistance during vibratory extraction will decline to roughly 25% of the values presented in Fig. 8.
Two simulations are performed to extract a pile from the soil.In the first simulation, the pile is pulled out without vibration, while in the second simulation, the vibration is applied.The simulation time of 7 min is chosen based on the field measurements of the pile vibratory extraction reported by van Dorp et al. (2021).The pile was extracted from the soil with a velocity of 40.5 mm/s, and the resulting pile shaft resistance is shown in Fig. 9.When the pile is 17 m embedded in the seabed saturated sand, the numerical value for the pile shaft resistance is 7.748 MN.This value is consistent with the results obtained from Equation ( 22) where the coefficient proposed by DNV is used.However, the present numerical study did not consider the soil consolidation effect, which could lead to higher shaft resistance.When the vibration with a frequency of 19.5 Hz is applied during pile extraction, the pile shaft resistance decreases to 1.8 MN from 7.748 MN, which is only 23 % of the initial shaft resistance.This finding is in agreement with that proposed value of the β method for dense sand (Jonker, 1987).

Table 4
Input parameters for sand (API, 2007;DNV, 2014;ISO19902, 2007).High frequency motion of pile induced by vibratory hammer, which is typically about 20 Hz, makes the surrounding soil particles to oscillate vertically and this causes a change in soil skeleton.As a result of the change in soil structure, the soil resistance is lost significantly, and a state of liquefaction is reached.This phenomenon happens when the acceleration of the soil particles (in the case of saturated sand) reaches a value of about 1.5 g (Viking, 2006).However, based on the graph obtained for the cyclic extraction of pile, it is not possible to judge whether the liquefaction has happened or not.Occurrence of liquefaction makes soil resistance approaching to zero.In fact, the pore water pressure in the sand increases, which reduces the effective stress and causes the sand particles to lose contact with each other.As a result, the sand loses its ability to resist deformation.In practice, it means the uplift force only should overcome the weight of the monopile and the vibro-hammer.However, it is visible in Fig. 9 (cyclic plot) that sand still shows some resistance after pile extraction has been triggered, suggesting that the soil has undergone degradation because of the cyclic loading.
To better understand how soil particles move and how pile extraction affects the soil, a comparison was made between pure pulling out and cyclic extraction of the pile.Fig. 10 shows the vectorial presentation of soil particle velocity in the vertical direction at different stages of extraction.The figure indicates that the soil particle velocity is positive (in the z-direction), which is due to the lifting of soil particles attached to the pile by friction while the gravity pulls them down, resulting in a velocity in the z-direction.Comparing the soil particle movements at the beginning of extraction, a larger region of the soil moves when the pile is pulled out without applying vibration, which can explain the higher required uplift force for this case.During the extraction process when the pile is partly extracted, the magnitude of velocity vector in the cyclic extraction is lower compared to that of the pure uplift case, as the soil particles experience both upward and downward velocities.During each cycle, pile experiences both pulling and pushing movement induced by the hammer.Therefore, the velocity of the soil particles surrounding the pile are influenced by this action and the particle velocity reduces in vertical direction.The region affected and undergoing movement is larger than the case pile is pulled out by crane uplift force, indicating that more particles in the horizontal direction are affected by the vibratory method.Even at the final stage of extraction, right before the pile is completely released, the affected region of the soil in the vibratory extraction method is large, but with a lower velocity magnitude in the vertical direction compared to the pile pulling out method.
The first and third stress components of soil in horizontal and vertical directions, respectively, are presented in Fig. 11 when the pile is extracted half of its embedded length.By comparing S11 for the two extraction methods, it can be observed that during pulling out the pile a region of soil around the pile experiences tensile stresses which initiates from definition of cohesion in MC model.The amount of this tensile stress is higher in shear pulling out the pile compared with vibratory extraction (with frequency of 19.5 Hz) implying that vibration reduces the strength of bounds between the soil particles.

Parametric study
In order to achieve a comprehensive understanding on the vibrohammer's features and the rate at which piles are extracted influence the process, a series of computational experiments is conducted using the basic specifications of the monopile mentioned earlier in Section 3.1.Table 5 summarizes the parameters used in this analysis.

Frequency of vibration
The vertical oscillation of the vibro-hammer is generated by counterrotating eccentric masses.The pile's ability to be driven into the soil is largely influenced by the resonance frequency of the vibrator-pile-soil system.This is due to the fact that the soil reaches the maximum vertical vibration velocity, and the pile penetration slows down significantly as the soil and pile vibrate in phase (Massarsch et al., 2017).At resonance frequency, the relative movement between pile and soil is small which results in an almost static pile-soil interaction (Massarsch et al., 2017).Liquefaction can be induced in saturated sand if the vibration frequency is high enough and the duration of vibration is sufficient.When extracting a pile using vibro-hammer, the pile is vibrated until it begins to move downward.Then, a crane applies an upward force to lift the pile and the vibrator together (Meijer, 2022).To gain a deeper understanding of how frequency affects the process of vibratory pile extraction, a variety of frequencies are numerically analyzed.The eccentric moment considered in this series of simulations is based on the hammer used in the field measurements (M e = 640 kg.m) and the reference pile is considered.Fig. 12 shows the pile shaft resistance during extraction.
Degradation of soil resistance which means reduction in pile shaft resistance is clearly observed for different frequencies considered in this study.This implies that in order to fully extract the monopile using the vibratory method, a partial upward force is required in addition to the force to overcome the weight of the vibro-hammer and the monopile itself.One could gain a more comprehensive understanding of the vibrohammer's effectiveness by taking into account the weight of the hammer In order to have better understanding of Fig. 12, a statistical analysis is carried out to acquire important features of each graph.These statistical values are presented in Table 6.The mean value of extraction with 19.5 Hz is lower than the other frequencies which indicates that lower uplift force is needed during extraction phase compared to applying other frequencies.The frequency utilized in field tests is close this value, although some fluctuations were applied during the extraction process due to the non-uniform nature of soil resistance along the depth.The standard deviation of the cyclic extraction with the frequency of 23.3 Hz is lower than that of other frequencies which means lower variations of resistance, i.e., the required extraction force is less spread around the mean.This suggests that even this frequency could be suitable during extraction phase.In order to provide a better insight on effect of frequency on pile shaft resistance, an ABAQUS filtration feature "Butterworth" is used.
A Butterworth filter is a type of signal processing filter used to remove noise or unwanted frequencies from a signal while preserving the essential characteristics of the signal.The Butterworth filter is designed to have a maximally flat frequency response in the passband, which means that the signal in the passband is not attenuated or distorted as much as possible.It achieves this characteristic by sacrificing the steepness of the roll-off in the stopband.The filter's performance is controlled mainly by the cutoff frequency.The cutoff frequency determines the frequency at which the filter starts to attenuate the signal.In ABAQUS, filtration is applied by setting cutoff factor.The cutoff factor refers to a dimensionless value that determines how fast the filter attenuates frequencies beyond the cutoff frequency.In Fig. 13, the plots of the pile shaft resistance versus the extraction time are filtered with a cutoff factor of 0.1.Implementing this approach has revealed that there is minimal disparity in outcomes when employing various vibration frequencies (i.e.,19.5 Hz,23.3 Hz,27 Hz,and 30 Hz).However, it is noticeable that the vibration frequency of 30 Hz can lead to a reduction in the pile shaft resistance during certain stages of the extraction process as compared to other three frequencies (see Fig. 14).

Eccentric moment
The new generations of vibratory hammers, such as the one employed for the retrieval of monopile foundations in decommissioning of the Lely wind farm, excel at generating high-frequency oscillations and can bear a consistent upward force during these oscillations.The cyclic force generated by the vibro-hammer is solely in the vertical direction as the moment created by the counter-rotating masses cancel out in the horizontal direction.The cyclic force can be calculated using the following equation: The centrifugal force of the vibro-hammer, F c , can be calculated using Equation ( 16).Selection of the vibro-hammer capable to produce enough centrifugal force is crucial in decommissioning operations to ensure that the extractability is reached, i.e., the hammer can trigger pile Table 7 presents a statistical analysis summary of the different centrifugal forces applied to the reference frequency of 19.5 Hz.Based on the mean and maximum values in the table, it can be concluded that using the CV-320-VLT-U hammer is not a suitable option for vibratory extraction of the reference monopile.The synchronized hammers can extract the reference monopile, with CV-640-VLT-U being the most  effective since it has the lowest mean value of pile shaft resistance during the extraction.The table results also suggest that using hammers with higher centrifugal force does not necessarily make pile extraction easier because a larger amount of centrifugal force results in higher oscillation amplitude.A higher oscillation amplitude means that a larger volume of soil attached to the pile moves during each cycle which requires overcoming higher soil resistance.To decrease the amount of centrifugal force, the vibration frequency can be reduced, but this may also result in soil degradation and liquefaction not occurring.Therefore, it is crucial to apply the appropriate amount of centrifugal force for effective monopile retrieval.

Extraction rate
The extraction rate was set to V = 40.5 mm/s based on the field measurements data (van Dorp et al., 2021).This rate has been varied to higher values to investigate its effect on pile shaft resistance in the present study.Fig. 15 shows the pile shaft resistance versus different extraction rates.It is observed that increasing the extraction rate decreases the efficiency of the vibro-hammer and requires larger force to overcome the soil resistance.This phenomenon occurs because the soil skeleton is not degraded due to the short duration of exposure to the cyclic loading.Therefore, finding the proper extraction rate appears to be important during the vibratory pile removal.
Table 8 shows the statistical values for different extraction rates.The statistical analysis of different extraction rates reveals that when the extraction rate is increased from 40.5 mm/s to 81 mm/s, the mean value of soil resistance increases four-fold, and the maximum resistance becomes 1.5 times larger.However, as the extraction rate is increased further, the increase in mean and maximum values becomes less severe.The table shows that the differences in mean and maximum values of the soil resistance between 3V and 4V extraction rates are smaller, suggesting that there is a limit to the effectiveness of vibratory extraction rate beyond which the higher rates do not produce satisfactory results.

Conclusion
The majority of decommissioned offshore monopile foundations have been cut either at the seabed level or below it, resulting in a significant amount of steel remaining permanently in the seabed.Total removal of monopile foundations is a more sustainable and environmentally friendly option, as it allows for the retrieved steel to be recycled and reduces carbon emissions.The application of vibratory extraction is a promising method.The limitations of using vibrohammers are relatively minor, especially after developing Vibro-Lifting-Tool (VLT) which eliminates the need for extra equipment to up-end or down-end the monopile before installation or after extraction, respectively.One limitation of the vibratory pile extraction lies in the capacity of the clamps.During the vibratory pile driving, the only force acting on the pile is its own weight and the hammer weight, and the installation is carried out with the pile in a state of free hanging.However, the crane operator regulates the penetration rate by adjusting the crane load to ensure the pile is installed vertically.In the extraction process, an uplift force is applied to the pile to remove it, emphasizing the significance of the clamps' capacity to maintain a secure grip.However, this limitation cannot restrict the application of vibrohammers.In the present study, extraction of offshore monopiles which is challenging as it involves consideration of highly nonlinear behavior of soil around the monopile is numerically analyzed.The built-in Mohr-Coulomb model in ABAQUS has some limitations in the analyses of large deformations, especially, when the soil is exposed to cyclic loading where soil particles experience both compression and tension and exhibit complex behavior.Despite the simplicity of the soil model, the   -There is minimal disparity in results of the soil degradation when various vibration frequencies investigated in this study (i.e.,19.5 Hz,23.3 Hz,27 Hz,and 30 Hz) are used during the vibratory pile extraction.However, it is noticeable that the vibration frequency of 30 Hz can lead to a reduction in the pile shaft resistance during certain stages of the extraction process as compared to other three frequencies, -Selection of vibro-hammer able to produce enough eccentric moment is crucial; otherwise, soil degradation will not happen, and the hammer fails to loosen the sand.On the other hand, selection of a vibro-hammer with much higher capacity than required functions inefficiently, because high eccentric moment which produces more centrifugal force results in larger oscillation amplitude.This in turn causes the pile to be displaced more during each cycle and more volume of soil is exposed to cyclic loading.Therefore, trying to degrade more soil volume is not reasonable by spending, or better to say, wasting more energy,    -The extraction rate plays a crucial role in soil degradation.A higher extraction rate means less time is given for the soil degradation, which may lead to a failure in reducing the soil resistance.
The MMC methodology suggested in this study is an effective approach for improving the precision of numerical modeling of pile extraction from dense saturated sand.This technique includes nonlinear features of stress-strain behavior and is relatively straightforward.Additionally, it is less computationally expensive compared to other available models that rely on soil constitutive models.The classical MC model with constant magnitudes of friction and dilation angles results in unrealistic deformations in soil that is modifiable by adopting the proposed MMC method.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
They used CEL technique to analyze the soil displacement around the casing during vibratory extraction and evaluated parameters like extraction speed, soil displacement and stresses.Experimental research on pile extraction using vibratory hammer is seldomly found in the literature.A series of onshore and offshore tests during Years 2016-2020 were conducted by van Dorp et al. (2021) to assess the

Fig. 6 .
Fig. 6.Total pile shaft resistance in dense sand for a pile 4 m in diameter.

Fig. 10 .
Fig. 10.Soil particle flow during pulling out and cyclic extraction (a) start of extraction, (b) half of the extraction and (c) full extraction.

Fig. 13 .
Fig. 13.Effect of different vibration frequencies on pile shaft resistance.

Table 1
MMC parameters adopted in the present study.

Table 5
Overview of the parameters.

Table 6
Statistical values for different extraction frequencies.
Note: the values of the statistical values are in KN.

Table 7
Statistical values for different centrifugal forces.
Note: the values of the statistical values are in KN.

Table 8
Statistical values for different extraction rates.
Note: the values of the statistical values are in KN.S.Salahshour et al.