Reliability-based load and resistance factor design model for energy piles

Energy piles have been popular globally with functions of both pile foundation and ground source heat pumps. Although several researches have been devoted to the probabilistic design and assessment of energy piles, the corresponding procedures are too complicated for engineers. As a simple variant of the reliability-based design method, the load and resistance factor design (LRFD) approach for the geotechnical design of energy piles is presented in this study. Firstly, the load-transfer model for energy piles is developed to investigate the effect of cyclic thermal loading on the pile settlement. Then, the LRFD procedures based on first-order reliability method and target reliability method are implemented into two different constrained nonlinear optimization problems, respectively. The proposed LRFD model for energy piles is demonstrated through an example pile and a series of parametric analyses.

Geothermal energy, as one of the renewable energies, is rich in China. To improve the energy utilization structure and protect the environment, the utilization of the geothermal energy has been promoted heavily. Ground source heat pumps (GSHP) system have been proved as a mature technology to use the shallow geothermal energy. In recent years, energy piles have been popular globally with functions of both pile foundation and GSHP. In practice, energy piles are typically precast or cast-in-place piles. The heat-transfer pipes are placed in the pile and made of high-density polyethylene. The diameter and wall thickness of the pipe are 10-40 mm and 2-4 mm, respectively. In cool climates or cold seasons, the heat present in the ground is typically extracted and transferred to the superstructure, while the heat is typically extracted from the superstructure and injected into the ground in warm climates or during hot seasons 1 .
Several researches have suggested that the thermal loading has a significant effect on the response of energy piles (e.g., [2][3][4], especially for their long-term performance (e.g., [5][6][7]. Thus, besides satisfying the design requirements for traditional piles, the thermal loading should be properly considered in the design of energy piles. Several analysis methods have been developed to investigate the thermomechanical behavior of energy piles, such as the empirical method (e.g., [8][9][10] ), load-transfer method (e.g., [11][12][13], and Full numerical method (e.g., [14][15][16]. However, it is still difficult to ascertain the pile performance in the geotechnical design of energy piles due to the uncertainties in soil parameters, applied load, design model, etc. In the traditional deterministic design approach, a conservative factor of safety (FS) is typically used to avoid the effect of these uncertainties, where those uncertainties are not explicitly quantified. For example, to prevent any potential thermally-induced damage to the building, the FS value of the in-situ energy pile was designed as 13 in Switzerland 4,17 .
As a more rational design approach, the reliability-based design method was recently adopted in the geotechnical design of energy piles [18][19][20][21][22] . In the context of reliability-based design, the performance of energy piles was investigated probabilistically. Based on the probabilistic method, the reliability index or probability of failure is typically calculated and selected as the measure of the uncertainties. By selecting a target reliability index or acceptable failure probability, the least-cost design that satisfies the safety requirement is typically selected as the final design. However, it is difficult to characterize the accurate statistics of these uncertainties due to the limited available data, which could further lead to an under-design or over-design of energy pile. To address those problems, a better model based on robust design method was proposed for energy piles, where the design robustness was considered to find a final design that is less sensitive to the uncertainties 23 . As a simple variant of the reliability-based design method, the load and resistance factor design (LRFD) approach is more popular in geotechnical practice (e.g., 24,25 ), whereas no relevant research has been reported on the energy piles. www.nature.com/scientificreports/ To this end, this paper aims to present a reliability-based LRFD method for the geotechnical design of energy piles. Firstly, based on the load-transfer method, a deterministic design model considering the long-term performance of energy piles is presented. Then, the LRFD procedures based on first-order reliability method and target reliability method are described, respectively, which are implemented as constrained nonlinear optimization problems in this study. By solving the optimization problems, the partial factors of each random variable can be evaluated. Focusing on the energy pile failure due to excessive settlement, the application of this LRFD model for energy piles is illustrated through an example pile. Furthermore, a series of parametric analyses are conducted to investigate the feasibility of the proposed LRFD approach.

Deterministic design model for energy piles
Load-transfer model. As a practical geotechnical design tool for traditional piles, the load-transfer method was firstly proposed by Coyle and Reese 26 to investigate their stress and strain distributions along pile. Recently, the load-transfer method was further adopted in the geotechnical design of energy piles 12 . In the load-transfer analyses of energy piles, the tri-linear curve (e.g., 12 ) and hyperbolic curve (e.g., 11 ) are typically selected as the load-transfer curves. Taking the Frank and Zhao's tri-linear curve 27 as the backbone curve, Fig. 1 shows the load-transfer curves that represent the relationship between side friction (t s ) [or base resistance (t b )] and pile-soil relative displacement (s). Note, the thermally-induced reloading and unloading paths of the load-transfer curves follow the Masing rules 28 . For fine-grained soils, the curve slopes can be determined as follows 27,29 : where K s and K b are the slopes of the first linear branches for the shaft and base load-transfer curves (MPa/m) [as illustrated in Fig. 1], respectively; E M is the Menard pressuremeter modulus (MPa), D is the pile diameter (m), and q s and q b are the ultimate side resistance and ultimate base resistance (MPa), respectively. It is noted that the values for q s and q b are typically calculated by the total stress method (α-method) or effective stress method (β-method) depending on the soil type.
In this study, the load-transfer model for energy piles is established using the matrix displacement method due to its good computational efficiency. As shown in Fig. 1, the pile is typically divided into N 0 equal elements. Each element and node are numbered from top to bottom. Similar to the traditional piles, the static equilibrium equation of energy piles is expressed as follows: where P(z) is the axial force at the depth of z (kN), T s (z) is the side friction at the depth of z (kPa). Considering the effect of temperature change, the compression deformation of pile section (ds) is given as: where {K S } e i is the stiffness matrix on the pile-soil interface of Element i, C is the pile's circumference (m), k s,i and k b are the secant slopes on the side and base load-transfer curves for a certain displacement (MPa/m), respectively, K h is the pile-superstructure interaction stiffness (MPa/m), and δ ij is the Kronecker delta (i.e., if i = j, δ ij = 1; otherwise δ ij = 0). Next, based on Eqs. (4)-(6), the external equilibrium is expressed as: where [K P ] is the total stiffness matrix of pile, [K S ] is the total stiffness matrix on the pile-soil interface, T is the column vector of pile node displacement, {P} = {P M , 0, . . . , 0, 0} T is the column vector of external loads, and {P T } = {EAα�T, 0, . . . , 0, −EAα�T} T is the column vector of thermal loads. Equation (7) can be solved using the incremental iteration method. It is noted the pile response at the end of previous loading is selected as the initial state in the calculation of next loading. For more details, readers can refer to some recent publications such as Fei et al. 11 , Luo and Hu 19,23 and Hu and Luo 18 .
Serviceability limit state design. Several researches have been conducted to investigate the role of thermal loading in the geotechnical design of energy piles. It was suggested that the effect of thermal loading must be considered in the analysis and design of energy piles 10 . Compared with the ultimate limit state, the performance degradation of the energy piles due to the piping need in the serviceability limit state was suggested to be more considered in the design (e.g., 6 ). This paper focuses on the excessive settlement failure due to cyclic thermal loading. The temperature cycle mode has been suggested in our previous study 18 . Taking the energy pile alternating summer and winter operation as an example, Fig. 2 illustrates a typical temperature cycle with a maximum temperature change of ± 10 °C. Thus, the limit state function of excessive settlement can be defined: www.nature.com/scientificreports/ where s lim is the limiting pile settlement (m), which is typically selected as 5% of pile diameter, s 1 is the pile settlement that can be computed with the developed load-transfer model (m).

Load and resistance factor design approach
Several methods are typically used for the back-calculation of partial factors in the load and resistance factor design (LRFD), including but not limited to first-order second-moment method (FOSM), Monte Carlo simulations (MCS), first-order reliability method (FORM) and target reliability method (TRM). Due to the merits of information on design point and parametric sensitivities 30 , both FORM and TRM are adopted herein to estimate the partial factors.
First-order reliability method. Low  where x is the vector of random variables (x i ), f is the failure domain, n is the unrotated equivalent standard normal vector of n i , R is the correlation matrix of the random variables, µ N i and σ N i are the equivalent normal mean value and standard deviation of x i , respectively, Φ[] is the cumulative distribution function (CDF) of a standard normal distribution, F(x i ) is the original nonnormal CDF evaluated at x i . As noted by Low and Tang 31 , β calculated by varying n is more efficient and robust than that determined by varying x. Hence, the n-space procedure 31 is adopted in this study. Figure 3a illustrates the FORM-based procedure to determine the design point (n*) for target reliability index (β T ). For a certain design scenario, the calculation of β is implemented as a constrained optimization problem 30 : where d is the design parameter, g(d, n) is the limit state function. Based on Eq. (11), a series of reliability analyses need to be performed by varying d until the calculated β is equal to β T . It is noted that the requirement that x ∈ f in Eq. (9) is imposed as a constraint of g(d, n) = 0 in Eq. (11).
Find: d Find: n Subject to: g(d, n) = 0 Objective: minimizing β Subject to: β = β T has been demonstrated to be efficient and reliable (e.g., [32][33][34]. In this study, the partial factors are also estimated using TRM. As shown in Fig. 3b, n* can be obtained by solving the following non-linear optimization problem: The 'fmincon' algorithm in Matlab is adopted herein to solve the aforementioned constrained nonlinear multivariable functions [i.e., Eqs. (11) and (12)]. Based on the solution results (i.e., n*) in Fig. 3, the partial factors corresponding to β T can be computed as follows: where x * i is the equivalent design point of x i , γ x i is the partial factor of x i , and µ x i is the mean value of x i .
Ethical approval. The manuscript has been prepared by the contribution of all authors, it is the original authors work, it has not been published before, it has been solely submitted to this journal, and if accepted, it will not be submitted to any other journal in any language.

Case study
To illustrate the proposed LRFD procedure for energy piles, an in-situ energy pile case 35,36 is adopted in this study for demonstration. The related parameters for pile and soil are listed in Table 1. The pressuremeter test data are used herein to estimate the ultimate side resistance and ultimate base resistance of pile 37,38 . As the thickness of silts is small, this study mainly focuses on the variability of sands. In addition, the pressuremeter test data is normalized with depth and the detailed statistics are further characterized. Thus, the mechanical load (P M ), pile-superstructure interaction stiffness (K h ), normalized pressuremeter modulus of sands ( E n M = E M z ), and normalized limit pressure of sands ( p n l = p l z ) are selected as the uncertain parameters. Table 2 shows the statistics of those random variables. The annual temperature change of the example pile follows the cyclic heating and cooling mode in Find: n Subject to: β = β T and g(d, n) ≥ 0 Objective: minimizing g(d, n) www.nature.com/scientificreports/ and critical parameters ( x = P M , K h , E n M , p n l ) on the long-term performance of energy piles. It is noted that the mean values of those random variables are used in the deterministic analysis. The number of pile element (N 0 ) in the finite difference model is set to be 3000 through a convergence analysis to ensure both the computational accuracy and efficiency. In practice, the pile settlement can be controlled by changing its geometric dimensions. Figure 4 shows the pile head settlements (s 1 ) at the end of 50 years along with the ratio of mechanical load to ultimate capacity (R m,u ) in different design scenarios. For energy piles, s 1 is induced by two types of loads, i.e., mechanical load and thermal load. It is worth noting that s 1 induced by the mechanical load is consistent with that of a normal pile with no temperature variation. It is obvious that both s 1 values decline with the increasing of L (or D) and finally stabilize. For example, the s 1 values with pile length of 9 m (R m,u = 92%), 12 m (R m,u = 55%) and 20 m (R m,u = 20%) are 18.0 mm, 5.7 mm and 4.3 mm, respectively. Similarly, when L = 12 m (Fig. 4b), the calculated s 1 values vary from 6.4 mm to 3.3 mm when D increases form 0.50 m (R m,u = 58%) to 1 m (R m,u = 21%). In addition, the s 1 values of a normal pile (i.e., black histogram in Fig. 4) are lower than those of an energy piles, especially when R m,u exceeds 60%. Hence, it can be easily concluded that the cyclic thermal loading has a negligible impact on the pile head settlement. The accumulated pile settlement due to cyclic thermal loading can be well controlled by ensuring that the mechanical load is under 60% of the ultimate bearing capacity of energy piles. Furthermore, considering 95% interval bound (i.e., two standard bounds) of P M , K h , E n M , and p n l (Table 2), the effect of those random parameters on s 1 is investigated. As shown in Fig. 5, s 1 values increase with the increasing of P M , whereas s 1 appears to be little affected by K h . In addition, larger E n M and p n l lead to a smaller s 1 .
Calibrating partial factors using first-order reliability method. Following the aforementioned FORM-based calibration procedure, a series of reliability analyses are firstly conducted. Considering the pile failure due to excessive settlement, the reliability index (β) at the end of 50 years is calculated for several combinations of L (10-14 m) and D (0.4-0.8 m). As shown in Fig. 6, β increases with the increasing of both L and D. For a given target reliability index (β T ), the most probable values of the uncertain parameters (n*) are available from the results of FORM-based reliability analysis. When β T = 2.33, the partial factors can be further calculated based on Eq. (14). As shown in Table 3, the calculated partial factors of mechanical load ( γ P M ), normalized pressuremeter modulus ( γ E n M ) and normalized limit pressure ( γ p n l ) vary in different design scenarios, whereas the change of the partial factor of the pile-superstructure interaction stiffness ( γ K h ) is negligible. The  www.nature.com/scientificreports/ γ K h value of 1.00 reveals that the uncertainty in K h has little effect on the pile failure due to excessive settlement, which is similar to the observation in the previous deterministic analysis (Fig. 5b).
Calibrating partial factors using target reliability method. As noted previously, in the FORM-based calibration procedure, d need to be varied by trial-and-error until the calculated β is equal to β T . In contrast, the target reliability method (TRM) is a more efficient method to back-calculate the partial factors because there is no need of varying d. Hence, based on Eqs. (12)- (14), the partial factors can be calculated by solving the constrained optimization problem [Eq. (12)]. It is noted that all the relevant parameters in the simulations are  www.nature.com/scientificreports/ taken from Tables 1 and 2. The effects of pile-superstructure interaction stiffness, operational mode of energy piles, and design parameter on the calibrated partial factors (i.e., γ P M , γ K h , γ E n M , and γ p n l ) are further investigated, respectively.
Following the cyclic heating and cooling load as shown in Fig. 2 (ΔT = 10 °C), a series of simulations based on the aforementioned TRM-based calibration procedure is performed. Figure 7 shows the calibrated partial factors for various target reliability index (β T ). It is obvious that a larger β T leads to a larger γ P M . This is mainly because a larger β T means a higher safety requirement. Moreover, a trade-off relationship between γ E n M and γ p n l can be observed with the increasing of β T , while γ K h remains approximately a fixed value of 1.00. As shown in Fig. 7, when β T = 2.33, the values of γ P M , γ K h , γ E n M , and γ p n l are 1.19, 1.00, 0.82, and 0.81, respectively. The calibration results using TRM are comparable with those using FORM ( Table 3).
Effect of pile-superstructure interaction stiffness on the calibrated partial factors. To further explore the effect of pile-superstructure interaction stiffness (K h ) on the calibrated partial factors, several other mean values [ µ(K h ) = 0.1 − 10 GPa/m ] and COVs ( δ(K h ) = 0.05 − 0.40 ) of K h are considered. Following the cyclic heating and cooling load as shown in Fig. 2 (ΔT = 10 °C), a series of simulations based on the aforementioned TRMbased calibration procedure is performed (β T = 2.33). As shown in Fig. 8, the change of partial factors with the increasing of µ(K h ) and δ(K h ) is negligible, which reveal that the partial factors of energy piles is less sensitive to the pile-superstructure interaction stiffness.
Effect of operational mode of energy piles on the calibrated partial factors. In the service life of energy piles, the pile performance is directly affected by its operational mode [operation time (t) and temperature amplitude (ΔT)]. Similarly, taking β T = 2.33 as an example and following the temperature cycle mode in Figs. 2, 9 shows the calibrated partial factors during 100-year operation time of the example energy pile. It is obvious that γ P M , γ E n M , and γ p n l change gradually with t, while the change of γ K h is negligible. In addition, when t > 10 year, the values of partial factors tend to stabilize, which reveal that the partial factors of energy piles is mainly affected by the thermal loading in the initial stage of service life.
To explore the effect of temperature change on the calibration results, several temperature amplitudes (ΔT ranges from 0 to ± 15 °C) are considered in the simulations. Similarly, based on Eqs. (12)- (14), Fig. 10 shows Table 3. Calibration results using first-order reliability method. L is the pile length, D is the pile diameter, β is the reliability index, γ P M is the partial factor of the mechanical load, γ K h is the partial factor of the pilesuperstructure interaction stiffness, γ E n M is the partial factor of the normalized pressuremeter modulus, γ p n l is the partial factor of normalized limit pressure. www.nature.com/scientificreports/ the calibrated partial factors at the end of 50 years for various temperature changes (β T = 2.33). It is obvious that the resulting partial factors are abruptly changed at the temperature change of ± 4 °C. The performance of the example energy pile is little affected by the thermal loading when ΔT is under ± 4 °C. As shown in Fig. 10, the values of γ P M , γ K h , γ E n M , and γ p n l when ΔT is under ± 4 °C are similar to these when ΔT = 0 °C (i.e., traditional pile). In addition, when ΔT is higher than ± 4 °C, the change of partial factors with the increasing of ΔT is negligible. It can be easily concluded that the calibrated partial factors are robust against the temperature change when the temperature amplitude of energy piles is is large enough (ΔT > ± 4 °C in this case study).

Design parameters Results
Calibrated partial factors in different design scenarios. The previous FORM-based calibration results (Table 3) revealed that the calibrated partial factors vary in different design scenarios. Hence, another series of pile lengths (L = 11-20 m) and pile diameters (D = 0.50-1.00 m) is considered in the simulations to explore the effect of design parameters on the calibrated partial factors. Taking β T = 2.33 and ΔT = ± 10 °C as an example, the calibrated partial factors at the end of 50 years can be obtained using TRM. As shown in Table 4, the changes of γ P M and γ K h with the increasing of L is negligible, while the values of γ E n M and γ p n l are abruptly changed when L > 13 m. Furthermore, it is obvious that the calibrated partial factors in each domain (i.e., L ≤ 13 m and L > 13 m) are similar to each other. Furthermore, the calibrated factors for various D are listed in Table 5. Similar to the observations in Table 4, the values of γ P M , γ K h , γ E n M , and γ p n l can be categorized into two parts: D ≤ 0.60 m and D > 0.60 m. It can be easily concluded that the design parameters have a significant impact in the LRFD of energy

Conclusions
This paper presents a reliability-based load and resistance factor design (LRFD) approach for energy piles considering the long-term performance degradation. Based on the load-transfer model, the serviceability limit state (settlement) of energy piles is investigated. The developed LRFD procedures are demonstrated to be effective through a case study. Based on the work presented, the following conclusions are drawn: (1) The LRFD model based on target reliability method (TRM) and first-order reliability method (FORM) are implemented into two different constrained nonlinear optimization problems, respectively. The partial factors calibrated using TRM are comparable with those using FORM, whereas TRM is an easier-to-use tool because there is no need of varying design parameters. (2) The operational mode of energy piles has a significant impact on the LRFD results. In the service life of energy piles, the calibrated partial factors change obviously in the initial stage and finally stabilize with time. Furthermore, the pile performance is little affected with the thermal loading when the temperature change is low (ΔT ≤ ± 4 °C), while the calibrated partial factors appear to be less sensitive to temperature amplitudes if ΔT is high enough (ΔT > ± 4 °C). (3) The partial factors of energy piles are robust against the pile-superstructure interaction stiffness in this case study. In addition, the calibrated partial factors vary in different design scenarios and the worst values of each partial factor are recommended as the final results in practice.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.