Exact Winkler Solution for Laterally Loaded Piles in Inhomogeneous Soil

A novel exact analytical solution is derived for the equation yð4Þ þ xny 1⁄4 0 in the region x ≥ 0, which is important for the analysis of piles in soil with stiffness varying with depth. To date, exact solutions for long piles are available only for the cases where n 1⁄4 −4, 0, and 1. For other values of the exponent n, solutions have formerly been obtained numerically, mainly by the finite-difference method or approximate analytical solutions. An inherent difficulty in obtaining solutions for long beams (which are used to model flexible piles) lies in the inability to isolate the regular, converging part of the solution over the singular part that diverges with increasing x. In this paper, an exact solution is derived for n > −4, focusing on the important case of semi-infinite beams. Key aspects of the problem such as the stiffness and flexibility matrices at the pile head, and the peak bending moments due to eccentrically acting lateral loads, are discussed. A novel approach for deriving Winkler spring moduli for combined force and moment loading is proposed and shown to provide good agreement with rigorous numerical continuum results. DOI: 10.1061/(ASCE)EM.1943-7889.0002125. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.


Introduction
The problem of a loaded beam supported by an inhomogeneous Winkler spring bed has been studied for over 100 years. The earliest analytical solutions known to the authors were provided by Hayashi (1921), who considered a linear inhomogeneity function for the variable Winkler modulus. The problem is still relevant in the modern day, where solutions can be employed to determine the deflection of a pile in an inhomogeneous soil under a lateral head load (Chang 1937), with solutions for power-law soil stiffness profiles being of particular utility. The problem is expressed by the equation y ð4Þ þ x n y ¼ 0, where x (≥ 0) and y are the independent and dependent variables, respectively, and n is a pertinent inhomogeneity parameter. Contrary to the corresponding problem of axial pile response in inhomogeneous soil which is straightforward to solve (Scott 1981;Guo 2012;Crispin et al. 2018), the problem of lateral pile response is harder to tackle. The underlying difficulties in obtaining analytical solutions lie in (1) the high order of the governing differential equation, (2) the variable coefficient (often associated with a noninteger power n), and (3) the divergent behavior of the solution at infinity. Hayashi (1921) provided two solutions for n ¼ 1. The first employed a Laplace transform to obtain an expression in an integral form. However, due to the complexities associated with evaluating (analytically or numerically) this integral, a second solution was provided as a formal power series. Similar power-series solutions have been reported by Rifaat (1935) and Hetenyi (1946); however, other than a similar solution attributed to Miche (1930) by Rifaat (1935), a workable integral solution was not provided prior to the seminal paper by Franklin and Scott (1979) using contour integration. Instead, numerical finite-difference schemes available in graphical form (Palmer and Thomson 1948;Gleser 1953;Matlock and Reese 1960) have been widely adopted in books and design manuals (Lam and Martin 1986;Reese and Van Impe 2011).
Despite the importance of the problem in engineering theory and practice and the intense research in the subject, until recently no significant progress in obtaining exact analytical solutions for n ≠ 1 has been achieved. The special case of n ¼ −4 has recently been solved by Froio and Rizzi (2016), who provided a constructive discussion on the physics of the solution. However, this is an inhomogeneous bed of limited practical importance for the pile problem considered here. Approximate analytical solutions obtained for arbitrary values of n derived by means of asymptotic methods and Wentzel-Kramers-Brillouin (WKB) perturbation series (Franklin and Scott 1979;Liang et al. 2014;Fradelos 2016) provide insight as to the attenuation of the solution with depth, yet do not offer dependable results for the stiffness and deflection at the pile head. Such results are available using the aforementioned numerical approaches or approximate analytical techniques based on energy principles and/or multilayer approximations of inhomogeneous soil (Mylonakis 1995;Mylonakis and Roumbas 2001;Mylonakis 2012, 2017). A summary of relevant energy methods and some refinements pertaining to dynamic loads is available in a recent review article by Mylonakis and Crispin (2021).
In recent years, interest in analytical treatment of the problem has re-emerged and some general analytical Winkler solutions encompassing arbitrary values of n have been formulated in terms of four generalized hypergeometric functions of the 0 F 3 type (Mylonakis 2004;Fradelos 2016;Crispin 2017;Froio and Rizzi 2017;Parashakis 2020;Rosendo and Albuquerque 2021;Cucuzza et al. 2022). These results extended the aforementioned classical power-series solutions, yet they are limited by their inability to handle the important case of long beams, which is used to model flexible piles. The problem lies in the divergent (singular) behavior of all four hypergeometric functions and the difficulty in extracting nondivergent (regular) solutions that tend to zero at infinity in the form of a linear combination of the four functions. Consequently, only piles of finite length can be treated, and specific boundary conditions at the pile tip should be considered. This limitation significantly complicates the solution and often leads to numerical instabilities for long piles (Fradelos 2016), providing the motivation for the herein reported work.
In the following, the problem is solved in an exact manner, with a focus on the important case of long piles (treated as semi-infinite beams) and arbitrary values of the inhomogeneity exponent n. To this end, the following steps have been undertaken: 1. The problem is formulated with the specific power-law inhomogeneity function defined and a novel definition of the Winkler parameter λ introduced. The desired form of the solution is detailed. 2. The existing solution for n ¼ 0 (homogeneous soil) is rederived to introduce the solution procedure used subsequently for long piles. 3. Two existing solutions for n ¼ 1 (linear soil profile) are compared and combined to develop a novel solution for long piles. 4. The power-series solution for arbitrary n (power-law soil profile) and short piles are rederived and extended to long piles through careful consideration of the asymptotic properties of the series. 5. Application of the Winkler solutions in Points 2, 3, and 4 to the pile problem is considered, including comparisons with numerical continuum solutions for the homogeneous and linear soil profiles. Regarding this final point, the main challenge in employing Winkler models lies in the assessment of the moduli of the relevant springs, kðzÞ. Various approaches are available to select kðzÞ, which is usually taken to be proportional to soil stiffness, E s (or G s ). Dependable expressions for kðzÞ from various sources have been summarized by Shadlou and Bhattacharya (2014), Anoyatis and Lemnitzer (2017), and Mylonakis and Crispin (2021). These are commonly obtained by calibrating the Winkler model against a more rigorous solution, usually a numerical continuum model implemented via finite or boundary elements. Calibrations of this type have been carried out by Biot (1937), Vesić (1961), Francis (1964), Yoshida and Yoshinaka (1972), Roesset (1980), Dobry et al. (1982), Gazetas and Dobry (1984), and Syngros (2004). Rigorous numerical continuum solutions for inhomogeneous soil are available (Banerjee and Davies 1978;Randolph 1981;Budhu and Davies 1988;Gazetas 1990;El-Marsafawi et al. 1992;Syngros 2004) that can be (and have been) used to calibrate the springs for inhomogeneous soils. However, the resulting Winkler modulus is dependent on both the soil stiffness profile as well as the pile head fixity conditions, thus requiring separate calibrations for different response parameters (i.e., head displacement, head rotation, and peak moment) and problem configurations (i.e., inhomogeneity parameters).
As an alternative, Mylonakis (2001) proposed a rational analytical approach based on a variant of the method by Vlasov and Leontiev (1966), originally developed for spread footings, that does not focus on calibrating a specific parameter. This solution was explored further by Karatzia and Mylonakis (2017), who approximated shape functions for inhomogeneous soil profiles using those determined from the Winkler solution when n ¼ 0. However, this was limited to only two boundary conditions: a pile under a head force without head rotation (fixed head) and a pile under a head force without head moment (free head); moment loading was not considered. In addition, the solution does not take advantage of the exact Winkler solutions presented here for inhomogeneous soils. In this work, the Vlasov and Leontiev (1966) solution is revisited to derive the Winkler spring coefficient for combined force and moment loading, properly accounting for the boundary conditions at the pile head. This solution employs the exact shape functions that are developed in the first part of this paper for inhomogeneous Winkler spring beds.
In summary, in this work, a solution is provided for the response of long piles in inhomogeneous soils under both lateral and moment head loading. This consists of (1) an exact solution to the relevant Winkler problem for the important case of semi-infinite beams, and (2) a systematic approach for obtaining the Winkler spring stiffness in inhomogeneous soil profiles for combined pile head loading. The proposed novel solutions can be used in engineering applications as well as for assessing other related models in geotechnics and engineering mechanics. Comparisons with alternative Winkler calibrations and numerical continuum solutions are provided.

The Problem
The classical Winkler model for a laterally loaded pile is shown in Fig. 1(a). The pile is treated as a cylindrical Euler-Bernoulli elastic solid beam of Young's modulus E p , length L, diameter D, and second moment of area I (¼ πD 4 =64). The soil reaction is modeled using distributed Winkler springs of stiffness kðzÞ, which varies with depth z and is measured in units of force per length squared (modulus of subgrade reaction).
In this paper, the following power-law Winkler stiffness profile is considered, shown in Fig. 1 where k D = Winkler spring stiffness at a depth of one diameter; n = power-law exponent; and z 0 = characteristic length corresponding to the elevation above the pile head that kðzÞ (if extrapolated above the ground surface) reaches zero. z 0 is shown graphically in Fig. 1(b) and is given by where k 0 = Winkler spring stiffness at the ground surface. This flexible formulation allows a variety of different soil profiles to be modeled. In particular, existing solutions are available for some special cases, including n ¼ 0, which produces a constant (homogeneous) stiffness profile (independent of z 0 , which is undefined for this case); n ¼ 1=2, which produces a profile with stiffness varying with the square root of depth, representative of some granular deposits (Darendeli 2001); and n ¼ 1, which produces a profile with stiffness varying linearly with depth often referred to as Gibson soil (Gibson 1974).
When designing a pile, of most interest to an engineer is usually the response at the pile head. Consider the combined lateral loading at the pile head shown in Fig. 1, consisting of a horizontal shear load P and a bending moment M. These actions result in a pile head deflection y 0 and rotation (anticlockwise positive) θ 0 , which can be related to the applied loading using either a (2 × 2) stiffness matrix K or a corresponding flexibility matrix F as shown in Eq. (2) where K 11 , K 22 , and K 12 (¼ K 21 ) = swaying, rocking, and crossswaying rocking stiffness coefficients, respectively; and F 11 , F 22 , and F 12 (¼ F 21 ) = corresponding flexibilities. These matrices are related by the simple expression Because both the pile and soil are modeled as linearly elastic, both matrices are symmetric by reciprocity (i.e., K 12 ¼ K 21 and F 12 ¼ F 21 ). In addition, the sign of each term depends on the overall sign convention chosen. For the sign convention shown in Fig. 1 (used throughout this paper), the off-diagonal stiffness terms K 12 and K 21 are positive, whereas the off-diagonal flexibility terms F 12 and F 21 are negative.
These matrices are obtained by solving the governing equation of the problem and applying the pertinent boundary conditions for each term. With reference to the problem in Fig. 1, the governing equation for pile deflection yðzÞ can be obtained from the wellknown strength-of-materials solution of an Euler-Bernoulli beam under the action of a distributed load, qðzÞ (Chang 1937) where the prime = differentiation with respect to depth z. For the sign convention in Fig. 1 (positive bending moments inducing negative curvatures) the relations among pile rotation θðzÞ, bending moment MðzÞ, and shear force PðzÞ are given by θðzÞ ¼ þy 0 ðzÞ ð 4aÞ MðzÞ ¼ −E p Iy 00 ðzÞ ð 4bÞ PðzÞ ¼ þE p Iy 000 ðzÞ ð 4cÞ Inputting the Winkler stiffness profile under consideration [Eq. (1a)] into Eq. (3) yields the following governing equation This can be simplified by introducing a Winkler parameter, λ (carrying units of length −1 ) such that y 0000 ðzÞ þ ðn þ 4Þλ nþ4 ½z 0 þ z n yðzÞ ¼ 0 ð5bÞ Eq. (5c) is a more general form of the commonly employed Winkler parameter for homogeneous soil (Hetenyi 1946;Franklin and Scott 1979;Scott 1981) (also discussed in the following section). This parameter can be interpreted as a characteristic inverse length ("wave number"), which enables significant insights into the nature of the problem and properly normalizes results. Note that the selection of D and k D as reference values in Eq. (5c) is an arbitrary choice; any other depth and corresponding stiffness would lead to the same λ. This is clear from the alternative form of Eq. (5c) provided in Eq. (5d), which uses k 0 and z 0 directly, neither of which depend on the choice of reference depth With the introduction of λ, only three dimensionless parameters are required to describe any problem modeled as in Fig. 1: λL, λz 0 , and n. This is evident from dimensional analysis, considering the five relevant problem parameters (k 0 , z 0 , n, EI, and L) and the two fundamental dimensions (force and length), yield three governing dimensionless groups (i.e., 5-2). The product λL incorporates both pile-soil relative stiffness and geometric pile slenderness and can be interpreted as a "mechanical pile length" (Di Laora et al. 2012), and the product λz 0 incorporates a pile-soil stiffness ratio dependent on k 0 , z 0 , and n and given by For zero surface stiffness, this term is zero, and only λL and n are required to describe the problem (as demonstrated subsequently in this work). Therefore, λz 0 can be interpreted as a dimensionless surface stiffness parameter.

Homogeneous Stiffness Profile
A reference application that deserves discussion is that of a homogeneous Winkler bed, i.e., when kðzÞ ¼ k. In this case the governing Eq. (5b) can be simplified to Eq. (6a) with the substitution n ¼ 0; the latter equation has the simple solution given by Eq. (6c) y 0000 ðzÞ þ 4λ 4 yðzÞ ¼ 0 ð6aÞ yðzÞ ¼ e −λz ½A 1 cosðλzÞ þ A 2 sinðλzÞ þ e þλz ½A 3 cosðλzÞ þ A 4 sinðλzÞ ð6cÞ where A m terms = integration constants to be determined from the boundary conditions at the pile head and base. Eq. (6b) is obtained from Eq. (5c) using the same substitution, which produces the wellknown Winkler parameter employed by many authors (e.g., Hetenyi 1946;Scott 1981;Pender 1993;Mylonakis 1995;Mylonakis and Gazetas 1999). The factor of 4 is normally included by convention to simplify the solution in Eq. (6c). Fig. 2 shows the pile head stiffness and flexibility as a function of dimensionless pile length λL for different base conditions (the plotted results for inhomogeneous stiffness profiles will be discussed subsequently). These are obtained using the pile head boundary conditions from the definition of each stiffness/flexibility term and three different pile base conditions: floating (i.e., zero load and moment), fixed (i.e., zero displacement and rotation), and hinged (i.e., zero displacement and moment). A systematic approach for obtaining these results is included in the Supplemental Materials and has been discussed by Crispin (2022).
From Fig. 2, it is evident that beyond a certain length, the boundary conditions at the pile base have negligible effect on the pile head response. It is therefore common to treat piles beyond this length, termed the active pile length L a , as long piles (i.e., semi-infinite beams). The first two terms in Eq. (5), associated with the negative exponent e −λz , approach zero as z increases, whereas the second two terms, associated with the positive exponent e þλz , are unbounded. Therefore, integration constants A 3 and A 4 must be set to zero for long piles, significantly simplifying the solution.
Various definitions for L a are available such as λL a ¼ 2 ffiffi ffi 2 p , suggested by Hetenyi (1946) for homogeneous soils (see also Di Laora and Rovithis 2015; Karatzia and Mylonakis 2017; Mylonakis and Crispin 2021 for a more recent discussion). However, regardless of definition, this category of long piles normally includes most piles installed in practice (with some notable exceptions such as large-diameter offshore monopiles) and is, therefore, the focus of this paper.
The pile head stiffness and/or flexibility terms can be determined from Eq. (5) using the appropriate pairs of boundary conditions at the pile head. For this purpose, it is convenient to introduce four basic response functions describing the pile displacement under certain conditions. In this paper, the following four functions are employed: • Y 1 ðzÞ, the displaced shape due to unit displacement at the pile head under zero rotation, • Y 2 ðzÞ, the displaced shape due to unit rotation at the pile head under zero displacement, • Y 3 ðzÞ, the displaced shape due to unit load at the pile head under zero moment, and • Y 4 ðzÞ, the displaced shape due to unit moment at the pile head under zero load. These functions can be obtained using any pair of independent solutions to the governing equation [Eq. (5)] that approach zero at large z, g 1 ðzÞ, and g 2 ðzÞ. These are herein referred to as the regular part of the solution. For homogeneous soil, these can be the functions multiplying the integration constants A 1 and A 2 in Eq. (6) [i.e., g 1 ðzÞ ¼ e −λz cosðλzÞ, g 2 ðzÞ ¼ e −λz sinðλzÞ], or any linear combination of them, and result in the expressions in Eq. (7) (Hetenyi 1946;Pender 1993). With the exception of Y 1 ðzÞ, these functions are not dimensionless or unitary. These are plotted in Fig. 3 (the plotted results for inhomogeneous stiffness profiles will be discussed subsequently) Based on these definitions and the sign convention chosen ( Fig. 1), the familiar stiffness and flexibility terms can be obtained directly (e.g., Pender 1993) as follows Once these have been used to establish the pile head deflection y 0 and rotation θ 0 , the overall displaced shape of the pile can be readily obtained from Eqs. (9a) or (9b) which can also be input into Eq. (4) to obtain the bending moment and shear force variation with depth.

Linearly Varying Stiffness Profile
Analysis of piles in other soil profiles requires repeating this process (notably solving the governing equation) for different kðzÞ profiles. To the authors' knowledge, the homogeneous case is the only commonly encountered soil profile that can be solved with elementary functions. However, solutions are available for soil stiffness varying linearly with depth, one by Hayashi (1921) [repeated by Rifaat (1935) and Hetenyi (1946)], who used an infinite powerseries solution, and another by Franklin and Scott (1979), who employed contour integration [also similar to an alternative solution by Hayashi (1921)]; both of these solutions are revisited here. Consider first a Winkler modulus function increasing linearly from kð0Þ ¼ 0, obtained from Eq. (1) by setting n ¼ 1 and z 0 ¼ 0. In this case, Eq. (5) can be rewritten y 0000 ðzÞ þ 5λ 5 zyðzÞ ¼ 0 ð10bÞ Hayashi (1921) and later Rifaat (1935) and Hetenyi (1946) proposed the following solution to Eq. (10): where the B m terms = coefficients to be determined from the boundary conditions; and the h m ðzÞ terms were given in the original work as a power series, which, following more recent solutions (Mylonakis 2004 where 0 F 3 ð; b 1 ; b 2 ; b 3 ; xÞ = generalized hypergeometric function of argument x and parameters b 1 , b 2 , and b 3 (Olver et al. 2010). The derivatives of these functions are provided in the Supplemental Materials for convenience. These functions are now (relative to when the solution was first proposed) much easier to compute, with efficient built-in functions available in many modern programming languages. For piles of finite length, the head stiffness and flexibility terms for different pile lengths and base conditions are plotted in Fig. 2, obtained in the same way as the solutions for the homogeneous stiffness profile [see Supplemental Materials and Crispin (2022)]. However, all four functions in Eq. (11) are unbounded at large argument, making them problematic for application to long piles. In principle, an arbitrarily large length, greater than the active length, can be employed to get suitable results; however, this approach is cumbersome and requires enforcing specific boundary conditions at the pile tip. Franklin and Scott (1979) introduced an alternative solution involving integrating a complex variable t along a contour c m , as shown in Eq. (12) [modified to use the definition of λ in Eq. (10c)]. Hayashi (1921) also proposed a similar solution derived using the Laplace transform.
Repeated differentiation of this equation yields the expression for the jth derivative Performing integration by parts on the fourth derivative as follows yields This shows that if c m is chosen such that the integrated part (first term on the right-hand side) vanishes, yðzÞ ¼ f m ðzÞ is a solution to Eq. (10). This naturally occurs when both limits tend to infinity with the same complex argument as a root of ð−1Þ 1=5 . The five contours with this property chosen by Franklin and Scott (1979) are shown in Fig. 4.
This results in five functions, corresponding to computing f m ðzÞ from Eq. (12) with each contour c m . These are easiest to calculate by taking the contours to follow exactly the limits of their defined region (dashed lines in Fig. 4). Each of these can be separated into two arms, the first approaching zero from one limit, and the second diverging from zero toward the second limit. Because the contours are drawn to use the opposite signs along each shared arm, the functions sum to zero at any z and are therefore dependent. Furthermore, from the symmetry shown in Fig. 4, the following relations between the functions can be derived in terms of each function's real, R½f m ðzÞ, and imaginary, I½f m ðzÞ, components: These properties have been discussed further by Franklin and Scott (1979). Franklin and Scott (1979) proposed the following general solution to Eq. (10) and showed that each of the four solutions are independent where the C m terms = coefficients to be determined from the boundary conditions.
Taking the contours to follow exactly the limits of their defined region, the jth derivative of f 0 ðzÞ and f 1 ðzÞ are given by [modified from Franklin and Scott (1979) Scott (1981) provided tabulated results for the real and imaginary components of these expressions obtained with a numerical integration technique (those results employed the dimensionless independent variable x ¼ 5 1=5 λz). Liang et al. (2014) proposed an approximate analytical approach for evaluating these expressions using a power series and the WKB approximation (Bender and Orszag 1978).
From the sign of the terms multiplying z, it is evident that f 0 ðzÞ approach zero as z increases, and f 1 ðzÞ is unbounded. In addition, f 0 ðzÞ and f 1 ðzÞ are uniquely defined as solutions to Eq. (10) and by four criteria on the function and its derivatives. Therefore, it is possible to obtain these functions as a linear combination of the generalized hypergeometric functions in Eq. (11) The power-series solution in Eq. (11) has the following convenient property at z ¼ 0: where y ðjÞ ðzÞ = jth derivative of yðzÞ with respect to z, as previously. Therefore, by setting Eq. (18) equal to Eq. (17) and letting j vary from 0 to 3, each of the B m terms and, consequently, analytical expressions for Eq. (16) can be obtained This combined solution can then be input into Eq. (15) and the response functions (as defined for the homogeneous case) obtained by enforcing the corresponding boundary conditions, i.e., employing the relationships shown in Eq. (7) with g 1 ðzÞ ¼ R½f 0 ðzÞ and g 2 ðzÞ ¼ I½f 0 ðzÞ, yields the four basic response functions for n ¼ 1 as follows where h 1 ðzÞ, h 2 ðzÞ, h 3 ðzÞ, and h 4 ðzÞ = hypergeometric functions in Eq. (11). It is worth stressing that although all four functions h m ðzÞ diverge with increasing z, the basic response functions in Eq. (20) tend to zero as z approaches infinity. These results complement and extend the solution of Franklin and Scott (1979). The preceding solutions are plotted in Fig. 3, and tabulated results for the functions and their derivatives are included in the Supplemental Materials. The stiffness and flexibility matrices can be obtained directly from the derivatives of these functions with the relationships in Eq. (8), leading to an explicit solution for the stiffness and flexibility matrices In which Each of the terms in Eq. (21) has the same dependence on λ as in Eq. (8) (for homogeneous soil). Yet, λ in this case is proportional to ðk=E p Þ 1=5 instead of an exponent of 1/4 for homogeneous soil [Eq. (10c)]. Therefore, for linearly varying spring stiffness with depth, the pile response has a higher dependence on E p and a lower dependence on k (and therefore E s ) than for a homogeneous soil. This trend naturally gets more pronounced with increasing n, as discussed subsequently in this paper.
For a nonzero surface stiffness k 0 , Eq. (10) must be modified as follows y 0000 ðzÞ þ 5λ 5 ½z 0 þ zyðzÞ ¼ 0 where z 0 , as introduced in Eq. (1b), is the elevation above ground level where kðzÞ (when extrapolated) would reach zero [ Fig. 1(b)]. This governing equation can be transformed to Eq. (10) with the simple substitution z → z þ z 0 . Therefore, the solutions in Eqs. (10) and (15)  Explicit closed-form expressions for these terms are possible; however, these are long and impractical and therefore not provided here. Naturally, for Winkler spring stiffness decreasing with depth, solutions for long piles are not available due to the stiffness becoming negative beyond a certain depth.
Consider the simplest case with zero surface stiffness (z 0 ¼ 0). Firstly, a simple monomial trial solution yðzÞ ¼ z s can be substituted into Eq. (5b) to provide the indicial equation From this result, it can be observed that for n greater than −4, a solution to Eq. (5b) in terms of Frobenius series of step (n þ 4) is possible when s ¼ 0, 1, 2, or 3 as follows where α m = coefficient of the mth term. Inputting this solution into Eq. (5b) yields the relationship between sequential α m coefficients as follows which can be used to derive the following four solutions to Eq. (5b) that can be used with Eq. (11a) The derivatives of these functions, obtained using standard differentiation rules for hypergeometric functions, are provided in the Supplemental Materials for convenience. When n ¼ 1, these expressions are the same as those in Eq. (11), and when n ¼ 0, they are solutions to Eq. (6). For piles of finite length, the head stiffness and flexibility terms for different pile lengths and base conditions are plotted in Fig. 2 for n ¼ 1=2 and n ¼ 1. These were obtained in the same way as the solutions for the homogeneous stiffness profiles [Supplemental Materials and Crispin (2022)]. The following important trends are worthy of note: • The K 11 and K 22 terms are strictly increasing functions of pile length for floating piles, and decreasing functions of depth for piles fixed or hinged at the tip. • A floating pile can never have higher stiffness K 11 and K 22 than a pile fixed or hinged at the tip. This, however, is possible for the cross term K 12 . • The opposite naturally holds for the flexibility terms F 11 , F 22 , and F 12 . • Beyond a certain pile length, the changes on pile head stiffness become asymptotic, and solutions for different conditions at the tip converge. Like in the case of linearly varying stiffness, the four hypergeometric functions in Eq. (26) diverge with increasing z. Accordingly, application of this general solution to long piles requires extracting the regular part of the solution, i.e., obtaining at least two combinations of B m terms such that Eq. (11a) approaches zero for large argument. This can be done by inspection based on the asymptotic properties of the hypergeometric functions to get the regular part of the solution, resulting in The basic response functions are obtained, as before, by substituting Eq. (27) into the relationships in Eq. (7) Tabulated results for these functions and their derivatives are included in the Supplemental Materials. The stiffness and flexibility matrices can be obtained in explicit form directly from these functions by using the relationships shown in Eq. (8) as follows where λ nþ4 ¼ k D ½ðnþ4ÞE p ID n [given in Eqs. (5c) or (5d) with k 0 ¼ 0]; ν ¼ 1=ðn þ 4Þ [given in Eq. (26e)]. The resulting head stiffness terms for different n values between 0 and 2, the range most likely to be encountered in practice, are provided in Table 1 with a precision of four significant digits. The preceding explicit solutions can be used instead of the (limited) tabulated numerical results in the literature to compute the pile head stiffness and flexibility for any power-law spring bed.
Nonzero surface stiffness can be incorporated, as previously, with the simple substitution z → z þ z 0 . The response functions must be recalculated from the relationships in Eq. (7) with g 1 ðzÞ and g 2 ðzÞ taken as g 1 ðz þ z 0 Þ and g 2 ðz þ z 0 Þ, respectively, from Eq. (27) [or any pair of response functions from Eq. (28) with the same substitution]. The resulting stiffness and flexibility terms are obtained, again, from the relationships in Eq. (8). For the important case of a square-root Winkler stiffness profile, where n ¼ 1=2, possessing finite stiffness at the surface, these are plotted in Fig. 6 alongside the result from the numerical method employed by Franklin and Scott (1979) [tabulated by Scott (1981)]. Evidently, the results of the analytical and numerical solutions are practically indistinguishable.
As the dimensionless surface stiffness k 0 and the associated dimensionless parameter λz 0 increase, it is expected that the pile head stiffness and flexibility terms approach those for homogeneous soil because the Winkler stiffness will vary less over the active pile length. However, this is not apparent in Figs. 5 and 6 because, despite the soil profile approaching being homogeneous for increasing z 0 , the λ value from Eq. (5c) used to normalize results does not reduce to the homogeneous value in Eq. (6b). Instead, it remains dependent on the n value initially employed, which no longer affects the soil profile. An alternative normalization approach must therefore be employed to observe the natural trend. One approach, shown in Fig. 7, is to use λ 0 to normalize the stiffness terms in the ordinates of the graphs. This stiffness is defined in Eq. (30) and is analogous to λ for the constant portion of the Winkler spring stiffness. Evidently, the pile stiffness and flexibility terms now approach those for homogenous soil with increasing surface stiffness

Assessment of Winkler Spring Stiffness
Application of the Winkler model to the pile problem for any soil inhomogeneity function requires calibration of the Winkler spring stiffness, kðzÞ. Therefore, the soil continuum must be considered in full, adding significant complexity to the problem. Existing solutions (discussed further in "Introduction" section) often employed boundary or finite-element analysis; however, in this work, an approximate analytical solution is employed. Considering where r = distance from the center of the pile; θ = aperture angle from the direction of loading; u and v = radial and tangential displacements of the soil, respectively; G s ðzÞ = soil shear modulus profile with depth; and η s = compressibility constant, which is a function of the soil Poisson's ratio, ν s . The compressibility constant η s depends on the choice of restraint employed to simplify the three-dimensional problem; in this case Mylonakis (2001) proposed assuming the additional vertical normal stress arising from pile movement, σ z , to be small, to get η 2 s ¼ ð2 − ν s Þ=ð1 − ν s Þ. Evidently, obtaining analytical solutions for the threedimensional problem described by the partial differential equation in Eq. (31) is far more complex than solving the governing equation of the simplified Winkler model in Eq. (5). Nevertheless, a rational simplification to Eq. (31), using information obtained from the solution to the Winkler model, can be employed to reduce the problem to two dimensions. To this end, the horizontal soil displacement vector is first expressed in separable form using the planar radial and tangential soil displacement components uðr; θÞ and vðr; θÞ (defined at an arbitrary depth), times a dimensionless unitary shape function of depth, YðzÞ=Yð0Þ, which is derived from the (nonunitary) basic response function YðzÞ in Eq. (28), upon division by its value at the pile head, Yð0Þ as follows uðr; z; θÞ ¼ uðr; θÞYðzÞ=Yð0Þ ð 32aÞ vðr; z; θÞ ¼ vðr; θÞYðzÞ=Yð0Þ ð 32bÞ Vlasov and Leontiev (1966) employed a similar approximation, but instead, they used a shape function in the direction perpendicular to the foundation-soil interface. Previous authors have used this approach for laterally loaded piles with a shape function for the attenuation in the radial coordinate, r (Vallabhan and Das 1988;Guo and Lee 2001;Basu and Salgado 2008;Basu et al. 2009). However, this approach results in a higher-order Winkler model of the Pasternak type with an extra term proportional to the second derivative of pile displacement. This can be interpreted either as an additional bed of distributed rotational springs, a membrane under tension connecting the Winkler springs at their base, or a shear beam connecting the base of the springs. This term adds significant complexity to solving the Winkler problem; to the best of the authors' knowledge, analytical solutions to the Pasternak model are only available for homogeneous soil.
By inputting Eq. (32) into Eq. (31) and integrating over the vertical coordinate, a reduced two-dimensional pair of equilibrium equations is obtained (Mylonakis 2001) where b is a constant (dimensions of length −1 ) encompassing the properties and response of the soil and pile medium with depth. For the common case of zero tractions at the soil surface and zero displacements at an infinite depth, b is obtained from (modified from Remarkably, Eq. (33) perfectly matches the well-known planestrain problem of Baranov and Novak (Novak 1974) governing the dynamic response of a horizontal soil slice under harmonic loading, with the b term retaining information from the suppressed third dimension being equivalent to a wave number normalizing the radial coordinate r or the pile radius D=2. The Baranov-Novak solution for the Winkler spring stiffness is (Novak 1974) where s = bD=2 (instead of the dimensionless frequency a o in the original model); q ¼ s=η s ; and K μ ðÞ is the modified Bessel function of the second kind and order μ.
Because b is dependent on YðzÞ which, in turn, depends on the Winkler parameter λ, the modulus k is present on both sides of Eq. (34) so, in principle, it must be determined iteratively. Karatzia and Mylonakis (2017) suggested that a good approximation can be obtained using only one iteration with an initial value of k ≈ E s . In addition, they provided a simplified version of Eq. (34) using the asymptotic forms of the Bessel functions for small argument pertaining to static conditions as follows where χ ¼ e γ =4 ≈ 0.445; η 2 s ¼ ð2 − ν s Þ=ð1 − ν s Þ; and γ (≈ 0.557) = Euler's constant.
The shape function used to calculate b, YðzÞ, depends on both the soil stiffness profile and the head loading conditions. For a fixed-head pile in homogeneous soil, using YðzÞ ¼ Y 1 ðzÞ from Eq. (7a) results in b ¼ λ ffiffiffiffiffiffiffi ffi 2=3 p ; similarly, for a free-head pile (under load only), using YðzÞ ¼ Y 3 ðzÞ from Eq. (7a) results in b ¼ λ ffiffi ffi 2 p (Karatzia and Mylonakis 2017). For a general combination of head load P and moment M; ðzÞ can be taken from the displaced shape given in Eq. (9b). For homogeneous soil, this can be resolved analytically resulting in Eq. (36), which is plotted in Fig. 8 The resulting Winkler spring stiffness k, obtained using Eq. (36), is also plotted in Fig. 8(a), normalized by the corresponding value for fixed-head piles, k F . Evidently, despite it being common to only use the fixed-head value (corresponding to an applied load to moment ratio of 2λ) (Fig. 8), there is significant variation in Winkler spring stiffness depending on the load to moment ratio. Fig. 8(b) shows the difference in a combined stiffness term ðP − λMÞ=y 0 when using the common assumption of the fixedhead Winkler spring stiffness compared with properly accounting for the load-moment ratio via Eq. (36). This combined stiffness was underestimated by around 25% for free-head piles with no head moment using the traditional approach and by as much as 60% for other load-moment ratios.
It is inconvenient to calculate k for each specific load-moment ratio to be considered, particularly for inhomogeneous soils, where analytical solutions to Eq. (33c) are hard to obtain. Instead, it is possible to take advantage of the property that only three terms are required to uniquely define both the stiffness and flexibility matrices in Eq. (2). Indeed, by establishing the spring moduli for a fixed-head pile (lateral force only, zero rotation), k F , an applied load only (free-head and zero moment), k P , and an applied moment only (free-head and zero lateral force), k M , the corresponding stiffness terms K 11 , F 11 , and F 12 can be obtained, and the remaining terms calculated as follows Here, k P , k M , and k F are determined using their corresponding b terms, namely b F , b P , and b M , which are, in turn, calculated from Eq. (33c) using functions Y 1 ðzÞ, Y 3 ðzÞ, and Y 4 ðzÞ, respectively, as the shape function YðzÞ. The b values from numerical evaluation of these integrals are given in Table 2 for homogeneous, square-root, and linear soil profiles, the latter two possessing zero surface stiffness. The resulting Winkler spring stiffnesses are plotted in Fig. 9 and compared against some existing solutions, including those of Francis (1964), who doubled the result of Vesic (1961) for a beam on an elastic foundation; Roesset (1980), who suggested the simple expression k ¼ 1.2E s ; and Syngros (2004), who provided separate expressions for k F and k P . In general, the iterative solution to Eq. (34) showed good agreement with the existing solutions, and the simplified approach proposed by Karatzia and Mylonakis (2017) in Eq. (35) provided very similar results. Also, for Eqs. (34) or (35), the differences in Figs. 9(a-c) are simply a stretch of the E p =E sD axis by a factor ðb=b F Þ ðnþ4Þ . Therefore, by applying this factor to  E p =E sD in existing formulas for fixed-head piles, approximate solutions for other load-moment combinations can be obtained. For nonzero surface stiffness (z 0 > 0), b values from numerical evaluation of the integral in Eq. (33) are shown in Fig. 10. By normalizing these results with λ 0 [ Fig. 10(b)] and letting the dimensionless surface stiffness increase, the results approached the values obtained for a homogeneous spring bed.
In summary, the approach detailed in this section to calculate kðzÞ consists of the following steps: 1. Determine Y 1 ðzÞ, Y 3 ðzÞ, and Y 4 ðzÞ from Eq. (28) using the approximation k ¼ E s . Use the λ definition in Eq. (5c) as well as the substitution z → z þ z 0 for nonzero surface stiffness. 2. Input each of these into Eq. (33c) as YðzÞ to get b F , b P , and b M , respectively. 3. Calculate k F , k P , and k M , using Eq. (35) and the corresponding b term. 4. Calculate K 11 , F 11 , and F 12 using the solution in Eq. (29) and the corresponding k term. Use the modified solution for nonzero surface stiffness. 5. Calculate the remaining stiffness and flexibility terms using Eq. (37). If high precision is required, the improved k values obtained at Step 3 can be input back into Step 1 (replacing k ¼ E s ) and the process repeated until they converge within a chosen tolerance. Eq. (34) can also be used in place of Eq. (35), although, as shown subsequently, neither of these changes are necessary to get good results. An illustrative example of applying this method to a freehead pile is presented in the Appendix.
The solution, shown here for power-law soil stiffness profiles, can be modified for other soil profiles (including multiple layers) by using appropriate shape functions in Step 1 and an appropriate Winkler solution in Step 4. The resulting kðzÞ=E s ðzÞ values are independent of depth, and this is a soil-structure interaction property depending on the full problem configuration.

Comparison with Continuum Results
The performance of the proposed solution for predicting pile head stiffness for homogeneous soil is shown in Fig. 11, where results are presented against available Winkler solutions and numerical continuum solutions (including some fitted functions). Employing Eq. (34) to predict k F , k P , and k M provided predictions of K 11 , F 11 , and F 12 in very good agreement with the continuum results as well as the recent Winkler spring stiffness calibrations for K 11 and F 11 by Syngros (2004).
The simplified one-iteration approach using Eq. (35) performed similarly, often indistinguishable from the results obtained from Eq. (34). The remaining head stiffness terms, obtained from Eq. (37) for both methods also showed very good agreement with the continuum results. Although the more rigorous solution from Francis (1964) provided better results for K 11 than the simplified finite-element based expression by Roesset (1980), both solutions highlighted the problem with using the fixed-head Winkler calibration for all head stiffness terms. This was most pronounced when calculating F 11 [ Fig. 11(d)], which has the highest dependence on Winkler spring stiffness, but is also evident in Figs. 11(b, e, and f).
The performance of the proposed solution for a linear soil stiffness profile [E s ðzÞ ¼ E sD z=D] is shown in Fig. 12. The results, employing Eq. (21), were similar to those for homogeneous soil. However, because each term has a lower power dependence on Winkler spring stiffness than the homogeneous case, the impact of differences in the Winkler spring stiffness employed was lower, and all results were more dependent on pile stiffness, resulting in less variation in both the Winkler and continuum results.
The predicted maximum moment, M max , for a free-head pile (M ¼ 0) and a fixed-head pile (θ 0 ¼ 0) are shown in Fig. 13. For the former, this was obtained from the Winkler model using Eqs. (4b) and (9b) with the corresponding shape functions for each soil profile and was negative in the sign convention shown in Fig. 1(a). For homogeneous soils, M max occurs at a depth z M ¼ π=4λ ≈ 0.79=λ, and for the linearly increasing soil stiffness [E s ðzÞ ¼ E sD z=D], M max occurs at a depth z M ≈ 0.96=λ (compared with numerical continuum results in Fig. 14). For the latter, the maximum moment occurs directly at the pile head and is the fixing moment, M F . This is obtained directly from the stiffness matrix using the expression M F ¼ PK 12 =K 11 .
The results for different Winkler spring stiffness calibrations are compared with functions fitted to available numerical continuum results in Figs. 13 and 14. In general, the Winkler model predicted higher peak moments than the continuum formulas for free-head piles, possibly due to the discretization in the continuum results missing the true peak. Results from this paper and fitted formulas by Syngros (2004) for both free-and fixed-head piles matched well. However, the calibration by Roesset (1980) (for fixed-head piles) overpredicted the maximum moments for the free-head case, particularly in the linear stiffness profile. The peak moment depth predictions for both profiles showed good agreement with the continuum results well, except for the function provided by Davies and Budhu (1986) for homogeneous soil, which was only fitted for E p =E s < 10 4 and diverged significantly above this value.

Conclusions
The Winkler model for the response of laterally loaded piled foundations embedded in an inhomogeneous soil with stiffness varying according to a power law with depth [Eq. (1)] has been investigated. To this end, the following conclusions can be drawn: • Existing solutions for spring stiffness profiles varying linearly with depth were extended to obtain solutions for semi-infinite piles. This included head stiffness and flexibility matrices [Eq. (21) and Fig. 5  conditions at the pile head, a novel approach was introduced where three separate k values were calculated to obtain the stiffness and flexibility matrices covering all cases [ Table 2, Fig. 9, and Eq. (37)]. This was shown to be a significant improvement compared with selecting a single k value (Fig. 8) and in much better agreement with numerical continuum results for predicting head stiffness terms and the peak bending moment .
Appendix. Illustrative Example Scott (1981) analyzed the results of a number of pile tests undertaken in the Arkansas River. For one of the tests, Pile 2, Scott (1981) proposed predicting the pile response using a linear soil stiffness profile. To demonstrate the method proposed in this work, this test will be used as an illustrative example for predicting the pile head displacement at a given load. The pile is a 16-m-long steel pipe with 0.41-m external diameter, 7.9-mm wall thickness, and a reported E p I of 69 MNm 2 . Based on knowledge of the site and the pile test results, Scott (1981) estimated the soil stiffness varies according to the following linear function of depth: E s ðzÞ ¼ E s1 z, where E s1 ¼ 35 MN=m 3 is the rate of increase in E s per meter depth. In addition, strain gauges showed that the bending moments were negligible at a depth of 4.6 m, indicating that the pile can be treated as a semi-infinite beam and the solution developed in this work can be employed. The load was applied at a height of 30 mm, so the applied head moment can be neglected, and only F 11 is required to predict the pile head displacement. This is obtained directly using k p , corresponding to inputting Y 3 ðzÞ from Eq. (28) into Eq. (33c). From Table 2, this results in b P ¼ 1.659λ, where λ is first estimated from Eq. (5c) [simplified for this case in Eq. (10c)] using kðzÞ ≈ E s ðzÞ to give λ ¼ 0.63 m −1 . Inputting these into Eq. (35) yields k P ðzÞ ¼ 1.9E s ðzÞ, which corresponds to λ ¼ 0.72 m −1 . A single iteration was deemed suitable here.
Finally, the result in Eq. (29b) [simplified for his case in Eq. (21b)] can be employed with these values to give F 11 ¼ 36 mm=MN. This means that at an applied load of 191 kN, the pile is predicted to deflect 6.9 mm. The prediction of Scott (1981) was based on a value of k ¼ E s , which is closer to values employed for fixed-head piles; therefore, the resulting predicted head displacement of 10 mm is not directly comparable with the value obtained here.

Data Availability Statement
All data, models, and code generated or used during the study appear in the published article.

Acknowledgments
The first author would like to thank the Engineering and Physical Sciences Research Council for their support (Grant No. EP/ N509619/1). The authors are grateful to the anonymous reviewers who provided constructive feedback that improved the work.

Notation
The following symbols are used in this paper: A m , B m , C m = integration constants determined from boundary conditions; a 0 (¼ ωD=V s ) = dimensionless frequency; b, b P , b M , b F = Winkler stiffness calibration constants; c m = complex contour; D = pile diameter; E p , E s , E sD = Young's modulus of the pile, (homogeneous) soil, and soil at one diameter depth; E s1 = rate of increase of soil modulus with depth; F = pile head flexibility matrix; n F m ðÞ = generalized hypergeometric function; F 11 , F 12 , F 21 , F 22 = pile head flexibility terms; f m ðzÞ, g m ðzÞ, h m ðzÞ = solutions to governing equation; G s ðzÞ = soil shear modulus; I = pile second moment of area; K = pile head stiffness matrix; K μ ðÞ = modified Bessel function of the second kind and of order μ; K 11 , K 12 , K 21 , K 22 = pile head stiffness terms; kðzÞ, k 0 , k D = Winkler stiffness variation and values the surface and one diameter depth; k P , k M , k F = Winkler stiffness piles under load-only, momentonly, and fixed-head conditions; L, L a = pile length and active length; MðzÞ, M = bending moment distribution and applied head moment; M max , M F = maximum bending moment and fixed-head moment; n = Winkler stiffness profile exponent; PðzÞ, P = pile shear force distribution and applied head load; t = complex variable; u, P = pile head displacement and load vector; u, v = radial and tangential soil displacement; YðzÞ, Y m ðzÞ = displacement shape functions; yðzÞ, y 0 = pile displacement variation and pile head displacement; z = depth below ground level; z M = depth of maximum bending moment; z 0 = elevation above ground surface such that kðzÞ is zero; ΓðÞ = gamma function; γ = Euler's constant; η s = compressibility constant; θðzÞ, θ 0 = pile rotation variation and pile head rotation, respectively; λ, λ 0 = Winkler wavelength parameters; ν = 1=ðn þ 4Þ; and ν s = soil Poisson's ratio.