Examination of the reduction in natural frequency of laterally loaded piles due to strain-dependence of soil shear modulus

Society ’ s transition to low-carbon energy sources has resulted in significant expansion in offshore wind tech- nology. Most offshore wind turbines (OWTs) are founded on large-diameter steel tubular piles driven into the seabed, termed monopiles. OWTs are subjected to large numbers of dynamic and cyclic environmental loads during their lifetime. To avoid fatigue or resonance issues, accurate characterisation of soil-structure interaction under operating and extreme conditions is paramount. There is a need to further improve understanding of the dynamic lateral response of piles after large lateral load cycling (storm events). This paper presents the results of lateral load field tests on piles with two different slenderness ratios, where the natural frequency was measured after a series of applied lateral load-unload events. 1D numerical models of the piles were developed and several operating parameters were varied to identify the parameters of importance in the nonlinear soil-pile interaction. Approaches to estimate the operating soil shear modulus both in-situ and post load-unload were trialled, and five subgrade reaction models were used to couple soil and pile properties, to ascertain the performance at predicting the measured frequencies from the experiments. Results suggest that predicted frequencies are highly sensitive to choice of model and degradation method used.


Introduction
Until recently, offshore engineering has been dominated by oil and gas infrastructure, which generally comprises large, heavy platforms founded on driven, slender piles. Vertical gravity loads dominate the design of these systems, and lateral loading due to wind and waves tends to be of secondary importance. The advent of offshore wind led to the development of larger, single piled systems known as monopiles and currently over 87% of installed wind turbines offshore are founded on these (Chortis et al., 2020;Wind Europe, 2018). Monopiles supporting OWTs experience significantly different loading characteristics to oil and gas jacket piles, with large lateral loads and overturning moments, and low vertical forces (Arany et al., 2017). The recently completed Pile-Soil Analysis (PISA) joint industry project Byrne et al., 2020aByrne et al., , 2020bMcAdam et al., 2020;Zdravković et al., 2020) has developed a new design approach for monopiles with a slenderness (length/diameter, L/D) ratio < 6, which is now widely adopted across industry. For piles with larger slenderness ratios (L/D > 6), the design of these still mostly relies on design procedures originally developed for the oil and gas industry.
The lateral response of piles is typically calculated using Winkler beam-spring models (Dutta and Roy, 2002;Prendergast and Gavin, 2016;Winkler, 1867), where the pile is modelled using linear-elastic beam elements (Clough and Penzien, 1993;Kwon and Bang, 2000), and the soil reaction is modelled using a series of non-linear decoupled springs. This approach is commonly referred to as the p-y method, where p denotes the lateral soil reaction and y denotes the lateral pile displacement at discrete locations along a pile. Application of Winkler-type models for estimating the lateral response of piles was first suggested by Reese and Matlock (1956). The most commonly adopted p-y approach for slender piles with L/D > 6 is formulated in the American Petroleum Institute (API) design code and has been used extensively in practice for the design of offshore structures. The API p-y curves for sand take the form of a hyperbolic tangent function, see Eq. (1), and were originally developed from a series of fourteen lateral load tests conducted on slender piles (O'Neill and Murchinson, 1983) under short-term monotonic (static) loading conditions.
where the ultimate soil reaction p u , at a given depth is dependent on the vertical effective stress, friction angle, and pile diameter; z is the depth below ground to the p-y curve in question; and k is the coefficient of subgrade reaction, a depth-independent parameter for a given soil density that controls the initial stiffness response of the system. The factor A accounts for static or cyclic loading in a rudimentary manner. Eq. (1) was developed based on a series of monotonic load tests, and cyclic loading is incorporated only by the use of the factor A. Offshore structures tend to experience loading that is predominately cyclic in nature, so several researchers have noted the importance of a more accurate characterisation of the response of these systems under this type of loading. There have been a number of field tests performed on laterally loaded piles to investigate cyclic loading, but these have primarily focused on fine-grained soils to date (Matlock, 1970;Reese et al., 1989). For sands, only a very limited number of cyclic field tests have been performed, and the API cyclic approach is largely based on a single case study at Mustang Island (Reese et al., 1974). The API approach adopts a value of A = 0.9 for cyclic loading, while for monotonic (static) loading the value of A decreases linearly with depth from 3 at mudline to 0.9 at ≈ 2.6D below ground level (D is pile diameter), below which it is constant. It should be noted that the use of the A term both within and as a multiplier of the hyperbolic tangent function in Eq. (1) results in both the static and cyclic API p-y curves having the same initial stiffness, but the cyclic curves have reduced ultimate resistance near the ground surface. Despite the limited empirical evidence and the lack of a rigorous theoretical basis, the API cyclic p-y method has been widely used to estimate the lateral response and natural frequency of offshore jacket piles as well as monopiles. Kallehave et al. (2012) demonstrated that the API approach tended to under-predict the soil stiffness for large diameter monopiles, resulting in a ≈ 5-7% underestimation of fundamental frequency of OWTs at the Walney offshore wind farm. With the rapid expansion of the offshore wind industry, there has been a need to develop new design methods specifically tailored for low-slenderness monopiles. The PISA joint industry project Byrne et al., 2020bByrne et al., , 2017aByrne et al., , 2017bDoherty et al., 2015;Zdravković et al., 2020) has led to the development of a new Winkler beam-based approach for monotonic lateral loading of monopiles, which offers significantly improved prediction of the lateral response of more rigid piles with L/D < 6. However, there remains substantial levels of uncertainty regarding the dynamic response of piles during load-unload-reload and post-cyclic loading.
One of the main issues with characterising soil-pile interaction lies with the strain-dependence and non-linearity of soil stiffness. In tandem with this, there exists significant uncertainty surrounding how to couple soil and pile properties in a comprehensive manner. Even in the case of small-strain elastic vibrations, characterising the initial stiffness response of piles is subject to errors. Prendergast and Gavin (2016) investigated this issue by conducting small-strain vibration tests on two piles and comparing the frequency results to models developed using a range of available subgrade reaction formulae. Five subgrade reaction expressions from the literature were compared, as developed by Biot (1937), Vesic (1961), Klopple and Glock (Elachachi et al., 2004;Okeagu and Abdel-Sayed, 1984;Sadrekarimi and Akbarzad, 2009), Meyerhof and Baike (Okeagu and Abdel-Sayed, 1984;Sadrekarimi and Akbarzad, 2009), and Selvadurai (Elachachi et al., 2004;Sadrekarimi and Akbarzad, 2009). Taking the small-strain shear modulus information from the test site at Blessington, Ireland, calculated directly by multi-channel analysis of surface waves (MASW) (Donohue et al., 2004), the results of beam-Winkler models incorporating soil-pile coupling stiffness from each of the five subgrade reaction formulae were compared with the experimental data from the two test piles. Results showed that the predicted frequency differed from the experimental results by between 16.6% and 31.5% for the first pile, and 3.9%-14.8% for the second pile, highlighting the significant disparity that exists in the small-strain stiffness estimation for pile-soil interaction alone. Due to the importance of accurate estimation of this operating stiffness, a finite-element model updating approach was subsequently developed by the same authors to identify the actual stiffness and mass of the soil contributing to the dynamic motion (Prendergast et al., 2019;Wu et al., 2018). Though useful, the question of how to accurately characterise the small-strain stiffness without the need for FE updating remains an open issue. Furthermore, these studies used the shear modulus profile estimated directly from MASW, so there is a question over how accurately this captures the operating stiffness for the tested piles where some overburden had been removed prior to testing.
For structures installed offshore, the harsh nature of the loading environment implies that piles will be frequently loaded beyond their small-strain elastic regions, and the resulting changes in system natural frequency and its temporal evolution are therefore of interest. This paper examines the effect of lateral load history on the system frequency of piles in order to understand the various mechanisms contributing to the problem. A well-known shear modulus degradation approach (Oztoprak and Bolton, 2013) is employed in this work to calculate the expected change in soil shear modulus due to applied lateral loading, in order to assess the impact of prior loading on the fundamental frequency. A field test on two piles subjected to lateral monotonic loading is undertaken at Blessington, Ireland, whereby the frequency of the piles is obtained before and after applying lateral loads of varying magnitudes. A range of approaches to estimate the in-situ and degraded soil shear modulus is trialled, and the results from applying the same five subgrade reaction methods are studied to ascertain their performance under the conditions investigated.

Experimental testing
The experiments in this paper were conducted at a dense sand test site in Blessington, Ireland, with properties analogous to offshore sites. The site has been the location for several full and pilot-scale foundation tests over the last number of years (Gavin et al., 2009;Gavin and Lehane, 2007;Igoe et al., 2011;Prendergast et al., 2015Prendergast et al., , 2013Prendergast and Gavin, 2016). The site comprises very dense, fine sand, and the relative density varies between 90% and 100%. The bulk unit weight of the sand is 19.8 kN/m 3. The sand is highly angular and has a constant volume friction angle of 37 • , and a peak friction angle varying between 54 • and 42 • as determined from triaxial tests (Tolooiyan and Gavin, 2011). The specific gravity of particles is 2.69 and the maximum and minimum void ratios are 0.73 and 0.37, respectively. The sand is partially saturated above the water table, which is located at approximately 13 m below ground level (BGL), and the degree of saturation above the water table is up to 75%. The D 50 of the sand grains varies between 0.1 mm and 0.15 mm. Ten Cone Penetration Tests (CPTs) were performed, which reveal an average CPT tip resistance q c profile ≈ 10 MPa at ground level to ≈ 17 MPa at 2 m BGL, increasing to ≈ 20 MPa at 7 m BGL.
Two 0.34 m diameter open-ended steel pipe piles, with a wall thickness of 14 mm, were driven to an embedment of 7 m using a 5 tonne Junttan hydraulic hammer. Each pile contained a mass of steel at its head. The sand around the first pile, P1, was then excavated until the pile was embedded 4.5 m, and around the second pile, P2, until the embedment was 3.1 m, to investigate the influence of slenderness ratio (L/D). While the diameter of the piles is much lower than those expected for offshore monopiles, the slenderness ratio is intended to represent piles that can be categorised as 'flexible' and 'stiff'. Pile P1 has a L/D of ≈13, representative of flexible piles as might be used for jacket foundations; and P2 has a L/D of ≈9, representative of stiffer piles, and is within range of existing installed monopiles in operation (Arany et al., 2017;Peder Hyldal Sørensen and Bo Ibsen, 2013). It should be noted, however, that many existing wind farms have monopiles with lower L/D ratios than those tested in this paper, and new deployments have significantly lower L/D (as low as 3 in some cases). The relative flexibility of the piles tested facilitates comparison of the reduction in shear modulus with depth along each pile, and so are useful in this context. A lateral loading rig was set-up at each pile to perform the lateral load tests, and each pile was fitted with accelerometers for dynamic testing. Fig. 1 shows a schematic of the pile arrangement for the lateral load and vibration testing on both piles (and numerical models, see section 3), and Fig. 2 shows a photograph of the site arrangement.
For the lateral load tests, each pile was loaded with the jack in increments that corresponded to approximately 10% of each pile's assumed ultimate lateral resistance. Each load stage was maintained for 5 min to allow for creep effects. The load (denoted as H in Fig. 1) was applied at a height of 0.4 m above the new ground line for each pile. Unload-reload loops were undertaken during the testing. The displacement and rotation of each pile was measured using Linear Variable Differential Transducers (LVDTs) and inclinometers, respectively, at multiple locations above the ground lines (see Fig. 2). The applied load was measured using a load cell. The lateral displacement and rotation of each pile as measured at each ground line is shown in Fig. 3 Prior to laterally loading each pile, a vibration test was conducted using the procedure described herein. Each pile was impacted a number of times using a modal hammer (PCB Piezotronics, 2018) and the output acceleration was measured by three accelerometers located along the pile shaft, as shown in Figs. 1 and 2. A Frequency Response Function (FRF) was developed for each impact, which is obtained as the ratio of the Fourier transform of the measured output to the Fourier transform of the input force time history. An example of the measured signals and resulting FRF from P1 is shown in Fig. 4, where Fig. 4(a) shows a segment of the force time-history, Fig. 4(b) shows a segment of the acceleration measured by the top accelerometer, and Fig. 4(c) shows the resulting FRF.
In addition to the dynamic tests undertaken prior to the application of lateral loads to each pile, subsequent dynamic tests were conducted at each load-unload stage, to characterise the dynamic response in the inelastic strain region. For P1, dynamic tests were conducted at two unload points corresponding to applied lateral loads of 150 kN and 270 kN, see Fig. 3(a) for both points. For P2, dynamic tests were conducted after unloading at 60 kN, 120 kN, 225 kN, and 300 kN, see Fig. 3(c). At each unload point, a number of impact vibration tests were performed, and the mean and standard deviation of each trial is shown for both piles in Table 1. In general, the frequency measured for each pile is observed to decrease after higher lateral loads have been applied, except in P2 at load stage 225 kN, which exhibits a higher frequency than at the lower applied load stage of 120 kN. The experimental data in this section will be used to appraise the performance of a number of approaches used to estimate the reduction in soil shear modulus with strain, and soil-pile interaction coupling in subsequent sections of this paper.

Methodology
In this section, the methods employed to develop numerical models of the piles and to estimate the in-situ, and reduction in, soil shear modulus G 0 are discussed. Two soil-pile interaction models are developed; (1) A non-linear Winkler beam p-y model for estimating the belowground lateral displacements of piles, which is used as an input to the G 0 reduction method; and (2) a dynamic Winkler beam model with linear soil springs to estimate the natural frequency of the piles based on the small-strain stiffness at each load-unload stage.

Model 1 -Non-linear lateral load (p-y) model
The lateral load field tests described in section 2 did not include any instrumentation below the ground line, therefore it is not directly possible to measure the pile deflection or rotation profiles BGL. In order to estimate the pile deflection at locations below-ground, a 1D Winkler  L.J. Prendergast and D. Igoe Ocean Engineering 258 (2022) 111614 beam p-y model was developed, as described herein. The pile is modelled using Timoshenko beam elements (Tedesco et al., 1999), discretised into 0.1 m segments. The elastic modulus of the beams is 210 GPa, the Poisson's ratio, v s is 0.3, and the Timoshenko shear factor is 0.5. The API static p-y model (API, 2007) was used to characterise the soil load-displacement behaviour, see Eq.
(1). The soil unit weight was assumed to be 20 kN/m 3 and the friction angle was estimated initially using CPT correlations (from the CPT tests undertaken at Blessington), and adjusted to provide a good match with the pile response at ground line. The friction angle profile used in the analyses is shown in Fig. 5, and a comparison between the ground line responses from the field tests and the 1D p-y analyses is shown in Fig. 6. It is important to note that the process for determining friction angle was solely for the purpose of matching the piles' mudline responses in order to estimate the below ground pile deflections, rather than for a rigorous analysis based on soil properties. The below ground pile deflections calculated using the 1D   p-y model for each load level, which are used in the subsequent G 0 degradation models, are shown in Fig. 7.

Model 2 -Small-strain dynamic pile-soil interaction model
The model developed to estimate the natural frequency of the piles is described herein. Each pile is modelled using the stiffness matrix method (Tedesco et al., 1999), where the piles are considered as Euler-Bernoulli beam elements supported on discrete, linear Winker springs (Dutta and Roy, 2002;Prendergast and Gavin, 2016;Winkler, 1867;Wu et al., 2018). A schematic of the numerical arrangement is shown in Fig. 1. Each beam element has four degrees-of-freedom, and the consistent mass and stiffness matrices used to model these are available in Kwon and Bang (2000).
P1 is modelled using 72 beam elements, each 0.1 m in length. These are supported on 45 spring elements, to model the 4.5 m of embedment. The soil within the pile is considered as an additional mass, by altering the effective cross-section of pile elements below the depth of the soil plug, taken as 5 m in length from the pile tip (as measured on site). The bulk unit weight of this internal soil is assumed as 20 kN/m 3 (Prendergast et al., 2019). This internal soil is not considered to add any additional stiffness to the system, so the second moment of area is unaltered. The mass of steel within the pile head is considered as a lumped mass with m t = 30.2 kg. P2 is modelled using the same information as for P1, except only 32 springs are used to model the contact soil, due to the lower embedded depth. The mass of steel within the top of P2 is calculated as m t = 18.2 kg.
The soil-pile interaction is considered by means of the coefficient of subgrade reaction approach (Ashford and Juirnarongrit, 2003;Prendergast and Gavin, 2016;Sadrekarimi and Akbarzad, 2009), whereby the pile and soil properties are coupled together, and stiffness coefficients are derived for each spring. There are a number of formulations available in the literature, each deriving from differences in the assumptions underlying their derivation. Two common approaches for pile-soil interaction were formulated by Biot (1937) and Vesic (1961), as shown in Eqs. (2) and (3). While normally used for static problems, various studies have applied these to small-strain dynamic modelling of pile-soil systems, e.g. (Ashford and Juirnarongrit, 2003;Prendergast et al., 2019;Prendergast and Gavin, 2016).
where E 0 is the small-strain Young's modulus of soil (kN/m 2 ), D is the pile diameter (m), EI is the flexural rigidity of the pile (kN m 2 ), and v s is the Poisson's ratio of soil. Biot's solution derives from the problem of an infinite beam with a concentrated load resting on a 3D elastic continuum, where a correlation between the continuum elastic theory and the Winkler model was obtained by equating maximum moments in the infinite beam. Vesic's solution was obtained via similar means, but in this case the maximum displacements (not moments) were equated to derive the solution.
In addition to the aforementioned models for pile-soil interaction, three similar models in the literature, which are normally applied to buried circular conduits such as pipelines, are trialled in the present work (Okeagu and Abdel-Sayed, 1984;Prendergast and Gavin, 2016;Fig. 5. Friction angle profiles used in 1D API p-y analysis. Once the selected subgrade reaction model is used to derive discrete Winkler spring coefficients, a numerical model of the pile soil system is developed by assembling the elemental mass and stiffness matrices into a global system, where M denotes the global mass matrix and K denotes the global stiffness matrix. The natural frequencies of the model can be derived by solving the Eigenproblem (Clough and Penzien, 1993), shown in Eq. (7). where

Estimating the in-situ small-strain shear modulus
The in-situ small strain shear modulus, G 0,in-situ , for the Blessington site was derived from MASW shear wave velocity profiles in tandem with using the CPT-based correlation by Schnaid and Yu (2007) as follows: where σ ′ v0 is the vertical effective stress, q c is the CPT cone tip resistance, and p a is the reference atmospheric pressure (taken as 100 kPa).
As the CPTs were undertaken at the original ground level and pile tests P1 and P2 were carried out after soil excavation to new ground levels, various approaches have been considered to account for the change in this ground level (and therefore stress level) on the G 0 adjacent to each pile. For the analyses conducted, the different approaches are termed G 01 -G 04 in the subsections below.

Method 1 -G 01
The test site at Blessington is heavily over-consolidated and very consistent in terms of relative density. It can therefore be assumed that the different elevations of the ground surface for each pile test should not significantly affect the soil response (i.e. a CPT performed at each level would provide approximately the same profile when plotted against depth BGL). For the G 01 case, it is thus assumed that the CPT q c profile is only dependent on depth below ground level (not elevation) and therefore each G 0 is calculated for each pile test using Eq. (8) assuming the CPT profile from the original ground level begins at the new pile ground level. Hardin and Black (1966) have shown that the small-strain shear modulus is dependent on the mean effective stress level, p ′ , such that:

Method 2 -G 02
where m is a material constant. Oztoprak and Bolton (2013) showed that the value of m in sands varies from ≈0.5 at small to medium shear strain levels, increasing to close to 1.0 at large strain levels. Based on this, the G 02 analysis case involves calculating G 0 at any given elevation from the CPT profile using Eq. (8), and applying a reduction factor based on the vertical effective stress as follows: where σ ′ v0,OGL is the vertical effective stress at that elevation assuming the original ground level. It is assumed that σ ′ v0 ≈ p ′ for heavily overconsolidated sand (i.e. the at-rest coefficient of lateral earth pressure K 0 ≈ 1).

Method 3 -G 03
A third analysis case for G 0 was considered where G 0 at a given elevation is calculated directly from the MASW profile (not corrected for effective stress level). This is arguably the most direct way to estimate G 0 and requires limited processing.

Method 4 -G 04
A final analysis case for G 0 is considered where G 0 at a given elevation is calculated directly from the MASW profile and corrected for the change in effective stress level using Eq. (10). Fig. 8 shows the derived MASW and CPT-based in-situ G 0 profiles using the various methods G 01 -G 04 plotted relative to the depth of embedment of P1 ( Fig. 8(a)) and P2 ( Fig. 8(b)).

Reduction in shear modulus under lateral loading
In order to account for the reduction in lateral stiffness of each pile due to prior lateral loading, a number of different models were developed based on the commonly used shear modulus degradation approach by Oztoprak and Bolton (2013). The approach relates the secant shear modulus (G) at a given shear strain (γ) to the small-strain shear modulus (G 0 ) using a hyperbolic equation as follows: where γ e is the threshold elastic shear strain, γ r is the reference shear strain at G/G 0 = 0.5, and a is a curvature parameter. Lower-bound, mean, and upper-bound values are provided by Oztoprak and Bolton (2013) in Table 2 based on a fit to a database of over 450 tests from the literature. For each node depth BGL, z, the average shear strain γ i(z) mobilised by the lateral movement of each pile y i(z) , is given by Matlock (1970) as shown in Eq. (12) where v s is the Poisson's ratio (assumed = 0.3). For all G 0 cases considered, the pile deflection y i(z) was calculated at each depth using the 1D p-y model described in section 3.1.1. In practice, the pile impact tests were undertaken after the pile had been unloaded from different load levels (H max,pre-test ). To apply Eqs. (11) and (12) using a fully loaded deflection profile (i.e. without considering any unload) would result in significant under-estimation of unload-reload stiffness. To incorporate this unload, two different approaches were compared; (i) linear unload, and (ii) Masing's type unload. Fig. 9 demonstrates both assumed mechanisms. The linear unload ( Fig. 9(a)) assumes the pile head unload curve follows a line with same slope as the initial stiffness. The Masing's unload ( Fig. 9(b)) assumes the unload curve follows that of the initial monotonic backbone curve, reversed and scaled by a factor of two (Muravskii, 2009). The procedure to calculate the operational G for these load-unload cycles is: (i) to apply the maximum load value for a given load stage in the p-y model to generate the pile head load-displacement curve, (ii) unload from this value using both the linear and Masing's approaches to calculate the plastic displacement, then (iii) to calculate what equivalent applied load Q app results in the same plastic displacement, which is subsequently applied in the p-y model to obtain the depth-dependant shear strains (Eq. (12)) and the degraded shear modulus (Eq. (11)). While this procedure is theoretically simplified, it has some desirable features, namely; a) Its simplicity allows for easy adoption within existing pile lateral analysis software. For example, the G 0 degradation can be calculated within spreadsheets and be applied using a stiffness modification factor (i.e. k-multiplier = G/G 0 ). b) It results in a depth-dependent stiffness reduction profile based on the pile displacement and soil strains (i.e. large reductions in stiffness near the ground surface where displacements are largest, and smaller reductions in stiffness at depth). c) It results in a strain-dependent reduction in stiffness where piles subjected to small loads within the soil's elastic range will not experience any degradation in stiffness, and piles subjected to large loads, which result in a highly non-linear response, will experience a significant reduction in stiffness.
For a more rigorous solution to the problem of post-cycling lateral response of piles, kinematic hardening elasto-plastic p-y models have been developed, which include important features relating to lateral pile behaviour such as gapping, ratcheting, and rate effects (see Beuckelaers (2017); Beuckelaers et al., 2017;Hededal and Klinkvort, 2010). However, due to their complexity and the difficulty in calibrating the required input parameters, these methods have not been widely adopted to date.

Analysis and results
The results of applying the modified G 0 profiles from sections 3.2.1-3.2.4 into the small-strain pile-soil interaction model from section 3.1.2 are presented in this section. The procedure is as follows: A given G 0 profile is estimated and used in conjunction with the degradation model in Eq. (11) to generate an operational G profile, using both a linear

Table 2
Parameters for hyperbolic degradation model (Oztoprak and Bolton, 2013 unload and Masing's unload approximation. This is subsequently inputted into each of the five subgrade reaction models and then inputted into the small-strain pile-soil interaction numerical model, to calculate the natural frequency of each pile after each load-unload stage. Three estimates are used for each degradation model curve; lowerbound, mean, and upper-bound, using the parameters as shown in Table 2. Each of the five subgrade models in section 3.1.2 are used to generate spring stiffness values for each pile model. The resulting natural frequency predictions (eigenvalues) from the pile-soil interaction model are compared to the experimental data from the tests undertaken on each of the two piles before lateral loading is applied, and after each load-unload phase. In order to ascertain which combination of G 0 profile estimate and subgrade model provides the closest prediction of the natural frequency of piles corresponding to the in-situ (pre-loaded) situation, it is necessary to compare the results obtained by the various G 0 profiles applied in conjunction with the subgrade reaction models. Fig. 10 shows the results of applying the four G 0 profile estimations for P1 with the five subgrade reaction models, prior to the application of any lateral loading. The 'preload' mean experimental frequencies from Table 1 are compared to the frequencies estimated using the small-strain pile-soil interaction model for each G 0 profile. From the results in Fig. 10, the G 04 profile (Fig. 10 (d)) combined with the Vesic subgrade reaction model provides the closest estimate to the experimental frequencies, with a percentage difference of 1.2%. This appears sensible as the Vesic model is directly applicable to infinite-beams, and has been used in modified form for pile-soil interaction analyses (Ashford and Juirnarongrit, 2003). It should be noted also that the K&G model combined with the G 01 profile provided a similarly close estimate to the measured frequency, as evident in Fig. 10(a).
The same analysis as shown in Fig. 10 is produced in Fig. 11, but this time for P2. The mean 'pre-load' experimental frequencies in Table 1 are compared to the predictions from each model using the applied G 0 profile and each of the subgrade reaction models. For P2, the Meyerhof and Baike (M&B) subgrade reaction model combined with the G 04 profile (Fig. 11(d)) provides the closest prediction to the experimental  L.J. Prendergast and D. Igoe frequencies, with a percentage difference of 0.4%. It is noteworthy that all remaining G 0 profiles result in an underestimation of the system frequency, shown by the column plots in Fig. 11(a-c) remaining below the experimental frequency line plot.
The previous analyses focussed on establishing which combination of G 0 profile and subgrade reaction model resulted in the closest prediction of the measured frequencies prior to the application of any lateral loading; with the result that the G 04 profile combined with the Vesic model performed best for P1, and the G 04 profile combined with the M&B model performed best for P2. The subsequent analyses focus on how these two resulting models for each pile perform at predicting the change in the natural frequency of both piles, when combined with the degradation model in Eq. (11), and two load-unload estimations (linear and Masing's-type).
In relation to the P1 tests, two load-unload stages at 150 kN and 270 kN were undertaken, with mean measured experimental frequencies of 17.5 Hz and 16.77 Hz, respectively (see Table 1). For both of these loadunload stages, three estimates of the G 04 degradation are obtained using Eq. (11), lower-bound, mean, and upper-bound; and two unload paths (linear and Masing's). These degraded G 04 profiles are shown in Fig. 12. Fig. 12(a) shows the in-situ and degraded G 04 profiles for the case of linear unloading at the 150 kN load stage. Fig. 12(b) shows the same data for the 270 kN load stage. Fig. 12(c) and (d) show the same data for the case of a Masing-type unload at the same load stages. The  degradation is depth-dependant in line with the deflection profile of the pile, and shows that P1 demonstrates flexible behaviour, due to the presence of two points of rotation.
The G 0 profiles in Fig. 12 are inputted into the pile-soil interaction model utilizing a Vesic subgrade reaction model and the results are shown in Fig. 13. Fig. 13(a) shows the results for the linear unload, and Fig. 13(b) shows those for the Masing's unload. The upper-bound degradation model provides the closest prediction in both cases. The pre-load experimental and numerical frequencies differ by 1.2%. When the 150 kN load was applied and released, the experimental and numerical frequencies predicted using the upper-bound degradation model differ by 3.5% for the linear unload, and 0.6% for the Masing's unload. For the 270 kN load case, the experimental and numerical frequencies differ by 21% for the linear unload, and 16.5% for the Masing's unload, suggesting the disparity grows with increased strain in the system. Both the linear and Masing's models over-predict the change in stiffness with strain, as demonstrated by the larger reduction in frequency than measured. This suggests that the method used may only be applicable to low to medium strain levels, or may require further calibration in the larger strain region.
In relation to P2, four load-unload stages at 60 kN, 120 kN, 225 kN, and 300 kN were undertaken, with mean measured experimental frequencies of 11.72 Hz, 10.82 Hz, 10.99 Hz, and 9.85 Hz (see Table 1). Three estimates of the degradation in G 04 (the best fit model to the preload frequencies), lower, mean, and upperbound are obtained for each load stage; and for two unload paths (linear and Masing's), as discussed for P1 above. The degraded profiles for the linear unload are shown in Fig. 14 (the Masing's unload are similar in shape but differ in magnitude); where Fig. 14(a) shows the data for the 60 kN stage, Fig. 14 The degraded G 0 profiles are inputted into the pile-soil interaction model using the Meyerhof and Baike subgrade reaction model, and the results are shown in Fig. 15. Fig. 15(a) shows the results for the linear unload, and Fig. 15(b) shows those for the Masing's unload.
Similarly to the case of P1, the upper-bound degradation model provides the closest prediction to the experimental frequencies across all load stages for the linear unload case (Fig. 15(a)), however the trend is not the same for the Masing's unload. For the pre-load case, the difference between experimental and numerical frequencies is 0.4%. For the 60 kN load stage, the upper-bound degradation model with the linear unload results in a 1% difference (over-prediction) between experimental and numerical frequencies, whereas for the Masing's unload, the average degradation model provides the closest estimate (0.3% difference). The same trend occurs at the next load stage, 120 kN, where the linear unload with upper-bound degradation results in a 1% difference in frequencies, and the Masing's unload with average degradation results in a 2.5% difference. For higher load stages, the trend re-aligns with that of P1, whereby the upper-bound degradation model with both linear and Masing's unload results in the closest predictions to the experimental frequencies in each case. Similarly to P1, the numerical predictions begin to deviate quite significantly from the experimentally measured frequencies at higher load stages, again suggesting that the overall reduction in shear modulus is over-estimated by the method used. There is additional uncertainty in the P2 case, however, as the measured frequency at the 225 kN load stage is higher than that at the previous load stage, suggesting some stiffness recovery may have occurred.
It should be noted that the higher load levels for P1 and P2 are approaching the pile lateral capacity and the soil is therefore strained far beyond the ranges for typical operational loads in practice. It should also be noted that in all of the analyses undertaken, the potential for added mass in the dynamic soil-pile interaction is neglected (Prendergast et al., 2019), however the results are still reasonable when compared to one another directly, as it is the relative difference that is of interest. Another important point to note is that the shear modulus degradation approach by Oztoprak and Bolton (2013), used in the analysis, is based on the secant stiffness of triaxial test data, rather than the unload-reload stiffness from cyclic triaxial tests, which could be deemed more appropriate for analysis of this problem. Nevertheless, the Masing's rule unload method with average degradation parameters provides a good estimate of the reduction in stiffness of both test piles at low to medium strain levels in this paper.

Conclusion
In this paper, the influence of various operating parameters in lateral soil-pile interaction on the natural frequency of piles embedded in sand was investigated. Two field tests were conducted whereby piles with two different L/D ratios were subjected to a number of lateral load-unload cycles of increasing amplitude. The natural frequency of each pile was measured using dynamic impact testing before application of the loads, and after the application of load-unload loops. A number of models were developed to characterise the operating G 0 profiles contributing to the small-strain stiffness of the two piles from both in-situ measurements of CPT tip resistance, and from MASW. The estimated soil-structure interaction stiffness using these G 0 profiles were inputted into five subgrade reaction formulae to estimate the operating soil-pile interaction stiffness. The degradation in G 0 with soil shear strain was estimated using a well-known degradation approach, with the severity of the degradation varied for three cases; lower-bound, average, and upperbound. Two different methods were developed for capturing the unload-reload stiffness; (i) based on a linear unloading, and (ii) based on unloading following Masing's rules. The resulting predictions from each model were compared to the experimental frequency measurements from both piles. The Masing's rule approach using average shear modulus degradation provided a good match with the experimental frequencies for both piles at small to medium strain levels. All of the methods over-predicted the reduction in natural frequency due to loading at large strains approaching the pile capacity. Future work will focus on the application of the simple stiffness degradation approach to p-y models.

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.