Experimental and Numerical Investigation of the Load-Bearing Mechanisms of Piles Socketed in Soft Rocks

Rock-socket pile design predominantly depends on the shaft resistance to support the load at the serviceability state. However, due to limited understanding of the pile–rock interactions, the pile capacity is normally calculated using empirical correlations. In this study, the load-bearing mechanisms of rock-socketed piles were investigated through a miniaturised pile–load test setup in a soft synthetic rock. X-ray CT imaging and numerical discrete element modelling were used to investigate the micro-mechanics influencing the load-bearing mechanisms at the pile–rock interface. The numerical pile model was established based on suitable constitutive models capable of simulating the soft rock behaviour. The analysis of X-ray CT images at various displacements revealed three different interface mechanisms, namely sliding, local shearing and progressive shearing. The numerical model validated this observed micro-mechanics in the rock asperities through the evolution of damage and micro-cracks. Insights from the experimental and numerical results indicated that the height of the rock asperities significantly dictates the failure mode. Results also illustrated that the shaft load–displacement response primarily depends on the forces acting on the leading edges of the pile asperity. In particular, it was observed that the bottom leading edge carried a predominant portion of the shaft loads due to its connectivity with the rock at the base of the pile. Though negligible, the forces on the trailing edges provided valuable information on the contribution of residual shaft resistance by the debris at the interface. Moreover, the numerical studies revealed the different failure modes at the pile–rock interface. The discussions presented in this study provide novel insights into the load-bearing mechanisms of piles socketed in soft rocks, which will help to improve design guidelines in the future. Reports the results of laboratory model-pile load tests incorporating X-ray CT imaging to study the load-bearing mechanisms of piles socketed in soft rocks. Experimentally investigates the roughness-dependent pile-rock interface behaviour, namely sliding, local shearing and progressive shearing with the help of 3D image visualization. Demonstrates the necessity of a user-defined damage plasticity model as opposed to conventional bonded particle models in numerical modelling of soft rocks. Correlates the interface mechanics with shaft resistance development along the entire length of the pile. Evaluates the asperity-level failure modes in soft rocks using DEM. Studies the effect of debris at the interface on the residual shaft resistance. Reports the results of laboratory model-pile load tests incorporating X-ray CT imaging to study the load-bearing mechanisms of piles socketed in soft rocks. Experimentally investigates the roughness-dependent pile-rock interface behaviour, namely sliding, local shearing and progressive shearing with the help of 3D image visualization. Demonstrates the necessity of a user-defined damage plasticity model as opposed to conventional bonded particle models in numerical modelling of soft rocks. Correlates the interface mechanics with shaft resistance development along the entire length of the pile. Evaluates the asperity-level failure modes in soft rocks using DEM. Studies the effect of debris at the interface on the residual shaft resistance.


Introduction
Rock-socketed pile (RSP) is preferred as the foundation solution for heavily loaded structures (e.g. high-rise buildings, bridges) when the bedrock is found within a reasonable depth. The design of the RSP requires a careful evaluation of the pile-rock interactions to satisfy the ultimate limit state and serviceability requirements. In designing the RSPs, the ultimate shaft and base resistances are essential for an accurate estimation of the axial load capacity. The performance of the RSPs has been extensively studied since the 1960s as researchers have tried to address various influential design factors (e.g. shaft roughness, smear materials, embedment length) that affect the design and construction of these piles. Earlier studies suggest that, in the serviceability state, the RSP capacity primarily depends on the ultimate shaft resistance due to its quicker mobilisation as compared to the ultimate base resistance (Whitaker and Wilfrid 1966;Haberfield and Collingwood 2006). The ultimate shaft resistance of the RSP can be correlated empirically to the UCS of the rock (Rosenberg and Journeaux 1976;Horvath and Kenney 1979;Meigh and Wolski 1979;Williams and Pells 1981;Horvath et al. 1983;Rowe and Armitage 1987;Carter and Kulhawy 1988;Kulhawy and Phoon 1993;O'Neill et al. 1996). Based on these field studies, the several key factors that affect the shaft resistance of the RSPs were identified as the normal stress at the pile-rock interface, the socket wall roughness, and the cleanliness of the interface. In addition to field pile load tests, several laboratory studies were conducted to explain the load-bearing mechanisms of the RSPs through direct shear tests on rock-rock and rock-concrete interfaces (Johnston and Lam 1989;Haberfield and Johnston 1994;Kodikara and Johnston 1994;Indraratna and Haque 1997;Indraratna et al. 1998;Seidel and Haberfield 2002). The interface shear strength was then established based around the fundamental models of Patton (1966) and Ladanyi and Archambault (1970). Furthermore, the slope failure analogue of Sokolovski (1960) was adopted by Seidel (1993) and Seidel and Haberfield (1995) to determine the peak shear resistance of regular-triangular rock asperities. However, the two-dimensional interface shear tests and models have their own limitations, such as the absence of tangential stress that is present around the shaft, and end of the socket boundary condition that affects the movement of debris generated during shearing. This prompted research on small-scale 3D pile-load tests in the laboratory to evaluate the load-bearing mechanisms of the RSPs. Dai et al. (2017) studied the effect of socket roughness on the shaft resistance through load tests performed on instrumented model RSPs. The ultimate shaft resistance was found to substantially increase because of the socket roughness, but an upper level of socket roughness was observed beyond which a minimal additional shaft resistance developed. Huang et al. (2018) conducted tests on model aluminium piles to study the effect of overburden pressure and socket length on the performance of the RSP. It was inferred that the range of the variations of ultimate shaft resistance under the influence of different rock socket length to diameter ratios was not more than 10%. The research also reported a very little to negligible influence of overburden pressure on the ultimate shaft capacity. Gutiérrez-Ch et al. (2021b) carried out centrifuge modelling of the model RSPs at 50 g to replicate geometric socket roughness and study the mechanisms at the shaft including the load transfer characteristics. However, the above-mentioned experimental studies were inadequate to capture the asperity-level mechanisms at the pile-rock interface (PRI) due to the absence of 3D visualisations. The first step in that regard was taken by Xu et al. (2020) who performed a miniaturised RSP experimental study and reported the mechanisms observed at the PRI subjected to incremental loading with the help of qualitative analysis of 3D X-ray CT imaging. Although this study made some progress to explain the asperity-level shearing mechanisms at the PRI, it could not fully evaluate the asperity-dependent (or roughness) failure modes. Therefore, it warrants an in-depth investigation of pile-rock interactions at the asperity-level including the 3D visualisation to understand the PRI shearing mechanisms for varying shaft roughness.
Apart from the small-scale laboratory tests, Discrete Element Method (DEM) has been widely utilised in recent years to investigate the behaviour of geomaterials at the grain-scale which provides valuable information about the grain kinematics, force and displacement fields, contact force network and fabric evolution (Shang et al. 2018;Saadat and Taheri 2020;Shen et al. 2020;Tran et al. 2020). The importance of grain-scale investigation in modelling the interface behaviour and the associated micro-mechanisms has been highlighted in several studies (Jiang et al. 2008(Jiang et al. , 2014O'Sullivan et al. 2008;Butlanska et al. 2014;Wang and Zhao 2014). Recently, Lambert and Coll (2014) indicated that it is possible to reproduce the experimental results of the direct shear testing of irregular rock interfaces using the commercially available DEM software PFC 3D . Furthermore, the modelling of a 3D rock-socketed pile was carried out by Gutiérrez-Ch et al. (2021a), where the authors utilised the material calibration parameters from the 2D direct shear tests to predict the shaft response. The authors proposed an empirical DEM-based factor to determine a lower bound value of the shaft resistance at displacements of about 1% of the pile diameter. However, capturing the PRI mechanisms without calibration with the field or laboratory pile load tests may affect their applicability in simulating the realistic boundary conditions. Also, their limited displacement criterion may not realistically represent the asperity shearing mechanisms as the shaft strength has not been fully mobilised. Therefore, they could not evaluate the debris movement at the PRI. Moreover, these research works focussed mainly on the ultimate shaft resistances and provide a very limited perspective on the asperitybased micro-mechanics at the shaft. In addition, the validity of the contact laws used in their DEM models was not verified against full pile-load test results.
Based on the above-mentioned review, this study addresses the following issues to further advance the rock-socketed pile research incorporating the micromechanics at the pile-rock interface.
1. The effect of tangential stress field at the PRI needs to be included for a realistic investigation of the pile-rock interactions influencing the 3D bearing failure; 2. The 3D visualisation of the asperity-level damage evolution at the PRI with incremental pile head displacements is to be incorporated in explaining the shearing mechanisms of PRIs with varying degrees of roughness; 3. The effect of rock at the pile base on the movement of debris generated from the shearing of asperities at large pile head displacements needs to be investigated; and 4. The influence of different DEM contact models and the associated asperity-dependent failure modes for modelling the full load-displacement behaviour of RSPs are to be evaluated.
This paper presents results of a comprehensive experimental investigation of model pile-load tests with X-ray CT imaging and 3D-DEM modelling to study the loadbearing mechanisms of the RSPs to address the abovementioned issues affecting the shaft resistance. Model piles with regular triangular asperities at the PRI were socketed into a soft synthetic rock and then subjected to axial loading with X-ray CT imaging. The pile-load test data for different asperity angles (0°, 7.5°, 17.5°, 22.5°) were analysed, and the micro-mechanics at the PRI were interpreted through qualitative analysis of the X-ray CT images. The peak shaft resistances were then evaluated with the existing empirical correlations. Finally, a numerical model of the RSP based on DEM modelling was established using different contact laws to study the micro-mechanics associated with the entire pile-load test response. This is the first research of its kind where the load-displacement response and associated asperitylevel failure modes are evaluated by employing model pile tests, X-ray CT imaging and DEM modelling.

Materials and Experimental Procedure
The experimental setup of the RSP was adopted from Xu et al. (2020), which consists of a miniaturised model pile, removable mould, base cover and pile test rig. Regular triangular asperities were used to represent the pile roughness in this study, and the dimensions of these piles were factored from the geometries used in the laboratory direct shear tests conducted on concrete-rock interfaces by Kodikara and Johnston (1994). They utilised a chord length of 54 mm to represent the 12.5° roughness profile which depicts the average chord length reported in the field (Pells et al. 1980;Seol et al. 2008). For this study, a fixed chord length of 9 mm (i.e. one-sixth of 54 mm) was used to simulate asperity angles of 0°, 7.5°, 17.5° and 22.5° around the socket wall (Fig. 1). These asperity profiles are also in line with the effective roughness heights proposed by Seidel and Collingwood (2001) through their correlations with the UCS of the rock. A 6 mm diameter aluminium pile was adopted to eliminate the influence of the radial boundary condition (Ovensen 1979) by limiting the size of the model.
In this study, the model piles were socketed into a soft synthetic rock which was produced from mixing of sand (0.45-0.6 mm), cement, kaolin and water in the proportion of 62.5:3.13:12.5:21.88 by weight. The mix proportion was established by trial and error to achieve a target soft rock strength and deformation properties (Table 1), which simulates a realistic field stiffness boundary condition (Pells et al. 1980). The normal stiffness boundary condition did not draw much attention in earlier small-scale pile-load tests (Dai et al. 2017;Xu et al. 2020). Therefore, the stiffness of the rock mass was kept constant and the synthetic soft rock properties were established based on material similitude analysis (Indraratna 1990) by maintaining the Poisson's ratio, critical strain and unconfined strength ratio between the field and lab approximately the same. This synthetic rock can be classified as an Argillaceous intermediate geomaterial as its properties fall on the boundary between clay and clay-based rocks (Hassan and O'Neill 1997). Since the stiffness of the rock mass has been transitioned from the field to the laboratory scale, the effects of rock mass jointing and quality have been ascertained in the synthetic rock samples.
The preparation of the model RSPs and experimental setup are illustrated in Fig. 2. The pile was first inserted into the base cover to check its rotatability and its resistance to swaying. The removable mould was placed over this setup, followed by the pouring of the rock slurry prepared using the HOBART® 5L mechanical mixer (Fig. 2a, b). The prepared samples were then completely sealed using a polythene wrap and kept in a standard curing chamber (20 °C and 95% relative humidity) for 28 days (Fig. 2c). After the curing process, the base cover was removed, and the RSP assembly was placed in the pile test rig for axial loading and X-ray CT imaging ( Fig. 2d-f). The load was applied on the top of the pile at 0.1 mm/min displacement rate using the INSTRON 5kN load frame. This displacement rate was selected based on tests conducted at different displacement rates (0.025, 0.05, 0.1 mm/min). Several factors, such as the X-ray CT imaging time (approx. 4 h/scan), age hardening and stress relaxation issues were taken into consideration to determine the optimum loading rate. For each roughness configuration, a minimum of three pile-load tests (two continuous and one staged) were conducted. The desired X-ray CT scan points were selected based on the evolution of the continuous load-displacement plot. The full image of the pile-rock system at different pile head displacements was taken using a high-resolution X-ray CT facility (ZEISS XRadia XRM-520Versa) available in the Department of Civil Engineering at Monash University. During the scan, the upward movement of the loading ram was locked. All scans were conducted at 46 microns voxel size using 140 kV and 10 W source setup with 7 s exposure time. To separate the shaft resistance from the total pile head load, a straight-cylindrical pile (0° asperity) was tested for the base resistance only, where the pile was rotated and removed from the socket, and then placed into the socket again to release the pile-rock bond and ensure negligible friction at the interface.

Load-Displacement Behaviour
The load-displacement responses for different pile shaft roughness subjected to vertical displacement are shown in Fig. 3a. It can be seen that the pile behaviour is initially linear elastic but exhibits strain hardening in the later stages. Also, the total load-carrying capacity of the pile increases with an increase in interface roughness. This total load-displacement behaviour encompasses the mechanisms happening at the shaft as well as at the base of the pile. Figure 3b illustrates the load-displacement responses of the 0 o pile with only the base bearing and a combination of the base and shaft bearings. It can be seen that the load increases nonlinearly with an increase in the vertical displacement of the pile. It is also evident that the shaft in this pile configuration supports almost no load, rendering the effect of the pile-rock bond and the frictional resistance at the interface negligible. Therefore, the load-displacement behaviour of the 0 o pile in this study can be considered as purely base resistance. This result, in turn, enables the determination of the shaft resistance for each asperity angle by subtracting the base resistance from the total load-displacement response (Fig. 4).
The results from Fig. 3a suggest that the experimental observations align with the findings of early researchers (Pells et al. 1980;Horvath et al. 1983) who reported a higher working load capacity for piles in sandstone and shale with  Experimental results on pile testing: a load-displacement behaviour with different roughness configurations, b influence of the shaft on the resistance of 0° roughness pile a higher roughness angle equivalent. This higher roughness attributes to an increment in the load-carrying capacity of the shaft due to the increased surface area between the rock and the pile at the PRI. As noted in Fig. 4, the ultimate shaft resistance of the pile increases by 176% when the asperity angle is increased from 7.5° to 17.5°, whereas the increase is only 20% when the roughness is increased from 17.5° to 22.5°. As the roughness of the pile is increased further, the rate of increase in pile shaft capacity declines and approaches an upper limit, beyond which it holds no significance. This phenomenon was also observed by Dai et al. (2017) through their experiments on model acrylic piles embedded into the rock. Therefore, the asperities, however tall after a particular limit, provide negligible additional resistance, as they are more susceptible to shear failure.

Failure Mechanisms at PRI
The PRI is susceptible to damage through the cylindrical expansion of the cavity accompanied by sliding and shearing of the asperities. The load-displacement behaviour of the pile moving through the rock socket is characterised by three different stages, namely: sliding (Stage I), local shearing (Stage II) and progressive shearing (Stage III). Figure 5 shows the total load-displacement behaviour of the 17.5° roughness pile for both continuous and staged loadings. The staged loading exhibits a comparatively stiffer response than the continuous loading due to the axial stress relaxation after each loading stage, as previously observed by Xu et al. (2020). For the same 17.5° roughness profile, the X-ray CT images acquired at each respective stage of loading are illustrated in Fig. 6. In the Stage I loading, the pile slides over the asperity face, which corresponds to the linear portion of the load-displacement curve (Fig. 5). This process imparts the load locally to the rock asperities through a squeezing mechanism. By inspecting the X-ray CT images in Fig. 6, this sliding can be visualised at the end of Stage I, where the voids start to open up radially at the PRI due to the gap created by the downward movement of the pile.
After sliding, the contact area between the pile and the rock at the PRI is restricted to a part of the asperity face, which causes a significant increase in the local contact stresses. Therefore, Stage II is characterised by localised crack formation and initial shearing of the rock asperities at the PRI, which gradually reduces the rate of increase of total load capacity, as seen from the change in slope of the shaft load-displacement curve ( Fig. 4). At this stage, the asperity shape is deformed to a maximum extent with yielding of the rock at the PRI.
Stage III of the pile-rock interaction mechanism comprises the progressive shearing of the compressed and locally sheared rock at the interface through the degradation of bonds in the shear plane of the asperity, and thereby a minimal increase in the socket resistance is detected. This stage releases debris into the socket, as evidenced in the XY section of the X-ray CT images (Fig. 6). The base resists a small portion of the total pile head load from the beginning, but when the asperities around the pile are completely damaged, the base starts to dominate the load-carrying capacity of the pile through a continuous strain-hardening approach. This observed behaviour of the rock contributing to the load-carrying capacity of the RSP correlates well with the mechanisms reported by previous researchers (Patton 1966;Ladanyi and Archambault 1970).
The shaft load-displacement response, as observed from Fig. 4, illustrates the mechanics at the PRI for different roughness configurations. In the behaviour of the 17.5° and 22.5° roughness profiles, it can be seen that no additional shaft response is available after the asperity has been predominantly sheared at the end of Stage III, as further process is characterised only by the movement of debris along the PRI. It can be evidenced by the almost constant load after the peak, which depicts frictional shearing. The subsequent strain hardening initiating after a vertical displacement of 90% of pile diameter can be attributed to the compression of debris in the socket, which results in an additional resistance at the interface. The rate of hardening after Stage III is higher in the 22.5° roughness pile due to the higher asperity height, resulting in a higher quantity of debris release into the socket. Similar mechanisms apply to higher roughness configurations, but the lower roughness profiles tend to behave differently due to a comparatively smaller asperity height. In that sense, the shaft response of the 7.5 o roughness configuration is considered here for analysis (Fig. 4). As characterised by sliding, Stage I constitutes the major portion of the entire load-displacement process, unlike the other roughness profiles. This phenomenon is due to the asperity height in this roughness configuration being insufficient to cause significant damage to the asperities. When the pile slides beyond 90% of the half asperity chord length, the socket resistance immediately starts to drop. The damage at the interface is minimal, and due to the rock stiffness, the compressed rock rebounds on the trailing edge of the pile asperities. This rebound pressure acts in the opposite direction to loading, thereby pushing the pile further into the socket. This constitutes the reason for the shaft capacity to reach the residual or lower limit at a vertical pile displacement of 8 mm. However, this phenomenon, as hypothesised from the X-ray CT images is investigated further in the final sections of this paper using a numerical approach.

Comparison with the Existing Design Methods
Many empirical correlations have been proposed in the past to determine the shaft resistance of RSPs based on field pile-load test results. The main form of these empirical relations can be generalised as where is the ultimate shaft resistance of the pile, and are the side resistance reduction factors which depend on rock strength and jointing, and q u is the UCS of rock mass. In this study, the empirical criteria that rely upon widely used UCS properties of the rock and adhere to the relationship in Eq. (1) are used for comparison with the experimental test results. For the evaluation purpose, the model piles with roughness configurations of 22.5° and 17.5° were classified as rough sockets of Class R4 and R3, respectively, and 7.5° and 0° as Class R2 and R1, respectively (Pells et al. 1980) as these asperity heights represent the scaled-down groove heights proposed for each roughness class by Pells  Table 2 shows the various shaft resistance formulations and their predictions of the ultimate shaft resistances for the soft synthetic rock used in this study. These values are compared with the experimental results in Fig. 7. Rosenberg and Journeaux (1976) utilised the results of Matich and Kozicki (1967) and compiled the field load test results in rocks with artificial roughening. The correlation provides a prediction between the upper and lower bound shaft resistances from this experimental study. It can be attributed to the highly fractured weak rock as reported in the literature, which should typically result in a lower shaft resistance as compared to the other higher roughness profiles (R3, R4). Meigh and Wolski (1979) reviewed the Horvath and Kenney (1979) relationships and predicted a lower bound value for the weak rock. This value is in good agreement with the results of the lower roughness configuration (R2). The same can be accredited to the results of Carter and Kulhawy (1988), who predicted the lower bound limits for the socket resistance in weak rocks based on the results of Rowe and Armitage (1984). Similarly, the formulations by Kaderabek and Reynolds (1981), Gupton andLogan (1984), Toh et al. (1989) and Stark et al. (2017) also yielded very conservative predictions. Therefore, these relationships are employed for predicting the lower bound value of shaft resistances. Rowe and Armitage (1987) carried out a comprehensive summary of the test results from over 20 field sites, and their predictions are the closest in line for the higher roughness configurations from our experimental tests (R3, R4) as the roughness classification in their study was also derived from Pells et al. (1980). A similar approach was taken by Zhang and Einstein (1998); however, their correlations over predict the shaft resistances due to their correlations primarily based on the various power-curve relationships and thereby affect the rocks in the lower UCS range. The experimental values for the higher roughness configurations (R3, R4) are in close agreement with the formulation by Williams and Pells (1981), who conducted tests in weak mudstone and supplemented the data from the literature for shale and sandstone. Even though the empirical correlation proposed by these authors is different from Eq. (1), the value of as 1 justifies its inclusion in these comparisons. Due to the usage of an aluminium pile in this experiment, the high mass modulus ratio yields the maximum value, thereby giving a higher shaft resistance despite the sockets having lower asperity heights as reported in Williams and Pells (1981). Horvath et al. (1983) proposed different relationships based on the roughness configuration of the rock socket. It closely predicts the shaft resistance value for the minimal roughness socket (R2). However, it underestimates the shaft resistances for the higher roughness configurations (R3, R4) due to the conservative predictions proposed in the study as the authors have reported the termination of some tests before mobilising the maximum shaft resistance. The empirical equations proposed by O'Neil and Reese (1999), Canadian Geotechnical Society (2006), and AASHTO (2008) ascertain the significance of discontinuities in rock mass but slightly overpredict the responses for very rough sockets (R4). The values are conservative for regular roughened sockets (R3) but in good agreement with the results of sockets with a lower bound value of roughness (R2). This discrepancy can be due to their formulations based on the recommendations of Horvath and Kenney (1979). Rezazadeh and Eslami (2017) compiled the empirical correlations of several field test results put forward by various researchers and proposed a general prediction irrespective of the shaft roughness. It is the reason for the predicted value to lie intermediate between the results of the almost smooth (R2) and the rough sockets (R3, R4). The predictions from Kulkarni and Dewaikar (2017) underestimate the shaft resistance of the roughened sockets and are more applicable for lower bound limits, but it has been reported by the author that the correlation is case sensitive and applies only to the site of study.
Correlations were also proposed by Seidel and Collingwood (2001), who considered the effect of construction parameters in the shaft resistance (Fig. 8). These correlations also contain the database of the field load tests by Kulhawy and Phoon (1993) for sockets of different roughness configurations. The experimental results from this study are consistent with the field responses compiled by Kulhawy and Phoon (1993) as they lie in the intermediate zone between the piles in rock and clay and thereby validate their property as an argillaceous intermediate geomaterial.
The experimental results are also contained in the predicted upper and lower bounds of the adhesion factors proposed by Seidel and Collingwood (2001) as the roughness values used in this study agree with the values obtained from the correlation of UCS and effective roughness height proposed by the same authors.
The base resistance of the pile in this experimental investigation is the load-displacement behaviour of the 0° roughness pile. As noted by Haberfield and Lochaden (2019), the minimum base resistance at a vertical displacement of 10% of the pile diameter should be a minimum of 5 times the UCS value of the rock (2.3 MPa in this case at a displacement of 0.6 mm). This value, when multiplied by the pile base surface area, gives a load of 64 N. It can be seen that this value is agreed upon precisely through the experimental results from this study (see Fig. 3b).
The experimental shaft and base resistance values from the miniaturised RSP tests tend to agree with the overall trend despite the variability of predictions from different  Table 2 with the experimental results from this study empirical formulations. It ascertains the applicability of the proposed methodology for studying the load-bearing mechanisms of the RSPs, which will be investigated numerically using the DEM approach below.

Methodology-DEM and Contact Model
Experiments with X-ray CT imaging on model RSPs are performed to study the mechanism at the PRI, but due to many influential parameters, the feasibility of conducting extensive experimental investigation is limited. Thereby the numerical models are adopted to calibrate the experimental response and then conduct parametric studies.
The DEM was first introduced by Cundall and Strack (1979) to capture the interaction of geomaterials through a system composed of rigid and finite-sized spheres. This system works on the principle of identifying the interaction forces between the particles based on force-displacement laws and then utilising Newton's second law to determine the accelerations of the individual bodies in the system. This process takes place for each timestep and is iteratively repeated until the end of the simulation. The interactions are created or destroyed by calculating the distance between two particles through the contact detection algorithm. The simulation progress is determined by the inputs given by the user in the form of a contact model which dictates the governing behaviour of each body inside the computational system.
In this study, simulations were carried out using a commercial DEM software PFC 3D (Itasca 2008). The software provides access to several inbuilt contact models, notably the Bonded Particle Model (BPM), to study the mechanical behaviour of cohesive materials (Potyondy and Cundall 2004). In the BPM, the behaviour of bonded interfaces is linear elastic till the bond strength is reached, and then it yields immediately by bond breakage (Fig. 9a), thereby making it computationally unsuitable for simulating soft rocks with a significant strain-softening response in the post-peak region. Several user-defined contact models exist for DEM to model the failure mechanisms of rocks (Šmilauer 2010;Scholtès and Donzé 2013;Onate et al. 2015). However, these models are unable to characterise the crack development due to mixed-mode failures. In this study, the cohesive damage plasticity model (CDPM) developed by Nguyen et al. (2017) was employed for simulating the mechanical response of soft rocks subjected to both normal and shear loadings (Fig. 9b). This model employs a damage law to account for the degradation of bonds between the particles, and thus the softening behaviour of normal and shear stresses can be obtained. In addition to this, a yield criterion regulating the yielding of bonds between the particles was developed based on both tensile and shear stresses. Hence, the model can handle both pure modes and mixed-mode failures. Also, the efficiency of this model has been upheld by earlier studies (Parisio et al. 2019;Saadat and Taheri 2019;Tran et al. 2021). Therefore, it serves as a suitable contact model to Fig. 8 Comparison of the experimental data with the field test results and correlations proposed and compiled by Seidel and Collingwood (2001) study the behaviour of piles socketed into soft rocks. The selection of the contact model as CDPM over the BPM for simulating the experimental response is explored further in the subsequent section.

Calibration of DEM Micro-parameters
The micro-parameters of the CDPM model for the soft rock need to be calibrated against its laboratory response to capture the physical behaviour. In this study, the UCS behaviour of the soft rock was utilised to determine the micro-parameters for the DEM simulations.
First, the selection of particle sizes needs to be carried out as it is an important aspect of the DEM simulation affecting the material characterisation (Potyondy 2015). A uniform particle size was adopted for the UCS simulations to account for the model homogeneity (Ding et al. 2014). The other micro-mechanical parameters were accounted for through the material contact parameters of the relevant contact model.
In the RSP simulation, the particles at the PRI should be sufficiently small and realistic to resonate well with the experimental behaviour. The particle size of 0.5 mm, similar to the D 50 of the sand in the experiment, was adopted at the PRI, thereby ensuring the effective interaction of the particles between the pile and the rock. However, simulating the entire model with this fine particle size would be inefficient, and therefore the particle scaling method was adopted to gradually increase the particle sizes radially away from the pile. This type of particle arrangement must ensure that the smaller particle must not migrate into the voids of the neighbouring regions. Therefore, particle sizes of 0.5, 0.75 and 1.125 mm were adopted by fixing an upscaling factor of 1.5. Unlike the BPM, it is imperative to perform separate UCS calibrations for the different particle sizes in CDPM as it has normal and shear softening micro-parameters that are dependent on the particle size. So, the particle size utilised in the UCS calibration must be similar for the pile simulations without any scaling procedure.
To calibrate the micro-parameters for the CDPM contact model employed in pile simulations, numerical samples for UCS tests were created by generating 15,682 balls corresponding to a constant porosity of 0.34. The micro-parameters (i.e. tensile strength, cohesion, effective modulus, stiffness ratio and softening parameters) were then identified by benchmarking the results of the UCS simulations against the experimental results, while the dilatancy and the friction angle were chosen to be 5.7° and 21.8°, respectively. The normal softening parameter was assumed to be equal to the shear softening parameter (Nguyen et al. 2017). A damping coefficient of 0.7 was selected based on the recommendations of Potyondy and Cundall (2004) to maintain equilibrium in the sample. In addition, to show the capability of the BPM model, a UCS test with a particle size of 1.125 mm was also numerically conducted for the same sample. Figure 10 shows the UCS behaviour obtained through the numerical simulations. It is evident that the DEM model with CDPM can capture the unconfined compression Fig. 9 Constitutive laws representing the failure modes: a Bonded particle model-BPM (Potyondy and Cundall 2004), b Cohesive damage plasticity model-CDPM (Nguyen et al. 2017) with the yield surface behaviour of a soft-rock specimen with good agreement to experimental results compared to the BPM model, which yields a highly brittle behaviour. The calibrated microparameters of the particles in the simulation, which correlate with the experimental results, are summarised in Table 3.

Numerical Simulation of the RSP
The numerical simulation of a large pile model requires massive computational effort and needs to be modelled appropriately to increase the feasibility. Earlier studies have tried to address the problem by adopting symmetrical sections to reduce the enormous calculations involved in the DEM analysis. For example, McDowell et al. (2012) and Falagush et al. (2015) simulated the cone penetration tests on sands in a 30 o model chamber and have reported no significant variations in the load-displacement behaviour when the system was downsized from 90° to 30° except for major undulations in the stress path. Working forward from this, Gutiérrez-Ch et al. (2021a) and Zhang and Fatahi (2021) have reported that segmenting the testing chamber of a pile-rock model into a 30 o and 45° section, respectively, increases the computational efficiency and preserves the quality of the model behaviour to a certain extent. However, no quantitative evidence has been reported to support this claim in rock media. To investigate this phenomenon, three different sized model sections (45°, 90°, 180°) were adopted to study the boundary effect, if any, arising from using narrow sections when the 3D pile behaviour is simulated.

Setup of the RSP Model
In the field behaviour of the RSPs, it is highly unlikely that the pile (represented by concrete) is susceptible to deformation during its interaction with the soft rock owing to the high concrete strength. Therefore, in this study, the pile was modelled by a combination of wall segments considering that it suffers zero deformation during loading. The discretization of the pile into several segments was carried out to analyse the load distribution across each segment with increased axial displacement. Figure 11 shows the sample preparation process of the RSP system in a 180° segment.
At first, frictionless walls were created to represent the pile and the outer boundaries of the vessel, as well as to support the generation of different particle assemblies (secondary walls) (Fig. 11a). Subsequently, three different assemblies of particles were generated in the spaces created by the pile, secondary walls and the vessel to present the rock mass, in which the particle size increased when radially moving away from the pile centre. This arrangement of particles is similar to the arrangement used in the study of McDowell et al. (2012) where cone penetration tests in soil were simulated in order to reduce the computational cost. The simulation was then run for 15,000 cycles to reduce the particle-particle overlaps. This process also reduced the movement of smaller particles into the adjacent regions due to the high magnitude of the locked-in stresses at this stage. Subsequently, the secondary frictionless walls were removed and then the model was allowed to reach the equilibrium stage. To get a better contact distribution and reduce these locked-in stresses, the radii of the particles were then  . Measurement regions were positioned inside the rock to calculate these isotropic stresses at each step during this process (Fig. 11b).
In the next step, the particles with less than three contacts were identified (floaters), and their radii were iteratively increased until all the particles in the assembly had a minimum of three contacts (Fig. 11c, d). Finally, the top wall was removed and the CDPM was applied to all the particle-particle contacts separated by less than 0.01 mm. The 'Linear model' proposed by Cundall and Strack (1979) was employed at wall-particle contacts and the effective modulus ( E * ) of these contacts was assigned a value twice as that of the contacts between particles of 0.5 mm size. The stiffness ratio of wall-particle contacts was assigned as 1.1. The wall-particle contacts representing the material boundary were assigned a micro-friction parameter of 0.4, similar to the particle-particle contacts in the simulation. It should be noted that contacts between the pile and particles had no friction, as the experimental results illustrated in Fig. 3b reveal that the friction parameter at the PRI is negligible.
After adding the gravitational load to the system, the loading process on the pile was simulated by vertically applying a downward velocity to the pile. This velocity was chosen as 0.02 m/s after performing sensitivity analysis and ensuring that the system was quasistatic.
The cross section and the results of the load-displacement test conducted on a 22.5° asperity pile for the three different sized model sections (45°, 90°, 180°) are shown in Fig. 12. It is evident from the result of the 45° segment that the boundary effect is a major concern while adopting narrow sections of the model to simulate the pile behaviour. When the contact area between the pile and the rock is increased, the difference in the load is minimised effectively by almost 30%. This result indicates that the overestimation of load resistance can be reduced significantly by adopting larger model segments.
It can be established from Fig. 12 that the response of the pile in the 90° segment closely matches with the 180° segment in the pre-peak region, but it shows a significant change in the strain-hardening response after the peak with an average difference in resistance of about 10%. It can be attributed to the boundary effect in the former, as the frictional boundary walls limit the movement of the particles in the lateral directions, thereby contributing to additional resistance in the post-peak region. To preserve the mechanism of the pile-rock behaviour and reduce the computational effort to some extent, the 180° segment was appropriately chosen to analyse the load-bearing mechanisms of the model pile.

Validation and Prediction of Pile-Load Tests
The impact of using CDPM over BPM in UCS tests was previously discussed in Sects. 3.1 and 3.2. A DEM sample was designed for the pile model using the 22.5° roughness configuration to evaluate the BPM and CDPM models. The  results of the pile simulations were compared with the experimental response, as illustrated in Fig. 13a.
It is evident that the BPM can almost reproduce the same linear elastic response just like the CDPM in Stage I (1 mm) but fails immediately after moving on to Stage II (1.6 mm) of the load-displacement behaviour. This is where the brittle nature of the BPM comes into effect; because when the contacts reach their yield point, their strength is abruptly reduced to zero. It results in the breakage of rock mass at the PRI, thus reducing the resistance to the shearing process. Therefore, Stage II of the load-displacement behaviour in BPM is characterised by a peak at 1 mm displacement, followed by a residual response that undergoes strain hardening due to the compression of debris in the socket. It finally undergoes tensile damage near the pile base which propagates diametrically, and the load resistance drops to a minimal value.
On the other hand, in CDPM, Stage II is characterised by pre-peak shearing, where the asperity undergoes damage, but the contacts do not yield completely. This benefit mainly arises from using the softening parameters in CDPM and validates its use over BPM in this study. It indicates that CDPM can characterise the macroscopic response of the RSPs in good agreement with the experimental results. This test served as the validation for the calibrated microparameters, and comparison with BPM further highlights the advantage of the CDPM model. Different configurations of the pile wall for the roughness angles of 0°, 7.5°, 17.5° and 22.5° were used to prepare the RSP samples (Fig. 13b). The results of the total load-displacement behaviour of the RSPs under different roughness profiles are shown in Fig. 14a.
It can be observed that the macroscopic responses of the DEM simulations exhibit a close agreement with the experimental results. The numerical simulation of the 22.5° roughness pile was first calibrated against its experimental response, and then the micro-parameters used in this calibration were employed to predict the behaviour of the other roughness configurations (0°, 7.5°, 17.5°). Figure 14b illustrates the shaft response of the RSPs with different roughness profiles for the numerical simulations and the experiments. It can be seen that the RSPs with higher roughness show a comparatively stiffer linear response in DEM when compared to the experimental behaviour at the very initial stages. This phenomenon is due to the near-perfect material preparation process in DEM and the immediate mobilisation of the socket resistance, unlike the experiment where the initial load is taken by the discrete cement bonds at the PRI and then by the asperity. The shaft resistance is mobilised quickly in the DEM simulation, whereas the experiment requires a certain threshold displacement to reach a stable pile-rock interaction process.

Asperity-Level Shaft Resistances of the RSP
To investigate the effect of the pile penetration into the rock socket, the pile in contact with the rock was divided into six different segments, each having half asperity chord length. Figure 15 illustrates the fractions of the total shaft resistance carried by the individual leading segments in terms of both the vertical and horizontal loads. The leading half segments (S2, S4, S6) are characterised as the ones establishing contact with the rock in the direction of pile loading (Fig. 15a). The remaining segments (S1, S3, S5), being on the trailing half of the asperity, do not experience any significant resistance from the rock during the downward vertical loading of the pile.
The normalisation of the loads with respect to overburden stress is negligible due to the small sample size involved in this study. The horizontal load development is observed to be higher than the vertical loads for all the roughness configurations. It is due to the inclined position of the asperity, which leads to the cylindrical expansion of the socket during its interaction with the rock.
It is interesting to note from Fig. 15a, b that the load carried by all the segments is almost the same till Stage I (1 mm) before the asperity deformation initiates. It is then marked by the increasing difference in the load distribution across the various pile segments. At the end of Stage II (1.6 mm), each segment starts following a different load-displacement path. The S2 and S4 segments experience a peak load similar to the total shaft behaviour and then exhibit a residual response. However, the S4 segment recovers quickly from the asperity damage, unlike the S2 segment and starts to contribute minimally to the overall shaft resistance. It should be noted that the S6 segment starts to exhibit a strainhardening behaviour and engages the major load from the shaft. This can also be observed from the behaviour of the S6 segment when the transition occurs from Stage I to Stage II. The cylindrical expansion of the socket during the displacement of the pile causes an increment in the horizontal normal stress on the S6 segment, thereby providing more resistance in the vertical direction.
The asperity shearing process is characterised by the formation of a shear plane arising as an impact of the high local stresses developed in the localised asperity region due to the squeezing of the rock at the PRI. The rock asperities are then degraded by the S2 and S4 segments and carried down into the shaft in terms of debris and therefore provide minimal resistance after the asperity damage. However, this is not the case for the asperity degraded by the S6 segment as it is connected to the base rock and does not have an opening for a shear plane to occur in the asperity. The inclined nature of the asperity at the base results in the transfer of the load imparted by the pile to transverse radially around the shaft and the base of the pile. Therefore, the residual shaft load-carrying capacity of the pile is governed significantly by the S6 segment. This leading half asperity segment also influences the behaviour of the shaft resistance for higher roughness configurations to adopt an almost elastic-perfectly plastic response (Fig. 4). Though the above discussion is valid for the S6 segment in all the roughness configurations, the behaviour of the segments S2 and S4 differ for piles with lower roughness (7.5°). By observing Fig. 15c, it is evident that the S2 and S4 segments exhibit a ductile response in the post-peak region as opposed to the brittle response in higher roughness profiles (17.5°, 22.5°). This variation is caused by the interface sliding having more influence than the asperity shearing on the pile-rock interaction.
As noted earlier, the shearing of asperity results in the expulsion of debris into the socket, which forms the basis of the residual resistance at the PRI. However, there is very minimal knowledge about the behaviour of debris during the vertical movement of the pile. To better understand this, the behaviour of the trailing edge segments (S1, S3 and S5) is analysed. These segments experience very minimal loads during the pile-rock interaction, but their load-displacement behaviour can give more insights into the development of the residual shaft response. It can be observed from Fig. 16a-c that the segments experience a compressive force in the opposite direction to loading. It is characterised by the compression of the debris in the rock socket. The S5 segment experiences more loads than S1 and S3 because the pile pushes the debris to accumulate near the base after its expulsion from the asperity face. This behaviour is observed for the 17.5° and 22.5° roughness piles (Fig. 16a, b). However, as the asperity height reduces, the amount of debris released also reduces. It results in the loads experienced by the trailing segments of the 17.5 o pile to be lower than the 22.5 o pile.
The debris deposit in the socket can be evidenced from the DEM sections for the 17.5° pile by observing the void formation at the PRI (Fig. 17). It can be seen that the white gap evolving at the PRI is effectively filled up by damaged particles at the end of Stage III (4.5 mm).
Regarding the 7.5° pile behaviour (Fig. 16c), it is interesting to note that all the trailing segments exhibit a similar resistance to the loads experienced during the pile-rock interaction. It is due to the very minimal release of debris into the rock cavity occurring because of a minor amount of asperity shearing. However, these loads are significantly higher in proportion to the total shaft load when compared to the higher roughness profiles. This is because the sliding of this pile over the asperity results in a significant elastic deformation of the cylindrical cavity, and when the asperity crosses the half chord length, the trailing face of the asperity comes into contact with the minimally damaged rock at the interface due to the elastic rebound of the rock, thereby experiencing highly proportionate loads. However, in the case of higher roughness profiles, the asperity deformation is mostly plastic.
It can be inferred from this analysis that the sliding and shearing mechanism at the asperity level plays a dominant role in influencing the load distribution across each segment of the pile. It should be noted that the load from the base is the same for all the roughness configurations, and it also mimics the load-displacement behaviour of the smooth pile (0° roughness).

Micro-Crack Development and Damage Evolution
The micro-crack formation in the rock was monitored during the application of vertical load on the pile, and the results are illustrated in Fig. 18. These micro-cracks occur due to the breakage of contact bonds between particles. It can be observed from Fig. 18a that the crack count progressively increases as the roughness of the pile increases. It is attributed to the additional shaft resistance experienced at the PRI due to increased asperity height. As the roughness parameter is incremented, the PRI tends to dilate more progressively but is constrained due to the stiffness of the surrounding rock. This dilation significantly increases the normal stresses at the interface, resulting in an upsurge of the micro-crack count. At the pile displacement of half asperity chord length (4.5 mm), this increase is about 65% and 246% when the height of the asperity increases by 0.445 mm (17.5°-22.5°) and by 0.82 mm (7.5°-17.5°), respectively. Also, it is interesting to note that the micro-crack count of the 7.5° and 0° roughness piles are almost the same. It directly relates to the interface behaviour, which evidently changes with asperity height because the PRI with lower roughness tends to fail mostly by sliding and, hence, exhibits minor damage at the PRI. However, increasing the roughness drives the asperity failure through damage of bonds, thus escalating the microcrack counts.
To support the notion further, the micro-crack count illustrating the various modes of failure of the contacts around the pile with respect to the pile displacement was investigated (Fig. 19).
The pile slides over the asperity till the cracks start to appear in the sample at the end of Stage I. It can be observed from the initiation of the micro-crack count from the zoomed section of Fig. 18a and the images of the DEM simulation in Fig. 18b. This behaviour at the PRI is then progressed by shearing and damage of the asperity, as evidenced by the change in slope of the curve at Stage II (Fig. 18a). The micro-crack count gradually increases at the PRI, and finally, most portions of the asperity are damaged completely at Stage III (Fig. 18b), where the unbonded particles occupy most of the shaft. The microcrack count is concentrated around the leading edge of the pile asperity as it is solely responsible for imparting the load to the rock. This process is characterised by the damage of the bonded contacts, which results in the formation of a shear plane through the rock asperity, splitting a significant portion of the asperity from the rock.
It can be illustrated from Fig. 19 that the shearing process evidently dominates the asperity failure mechanisms of the RSPs. This can also be visualised from the trend of the shear micro-crack count curve which resembles the results of the respective pile roughness configurations in Fig. 18a. Next to pure shear, a combination of shear and tensile failures, also known as the mixed mode, is responsible for the majority of the contact damage. It is also evident that the pure tensile failure is minimal irrespective of the pile roughness configuration. It should be noted that it is almost negligible for the 7.5° pile, and the asperity damage, however minimal, occurs mainly as a consequence of shearing. The evaluation of microcrack count ascertains the fact that asperity shearing is a predominant component of the failure mechanisms and dictates the development of interface resistance. The displacements associated with the micro-crack development may vary based on the asperity geometry.
On the other hand, at the base of the pile, the damage is characterised by huge compressive forces acting at every stage of pile displacement. These analyses depict their relevance to the mechanisms interpreted from the qualitative study of the X-ray CT images and hence support the experimental micro-mechanics.

Conclusions
The design of the piles socketed into rocks heavily depends upon the resistance along the shaft particularly in the serviceability state. In this study, an experimental investigation was conducted to study the load-displacement behaviour of the rock-socketed piles based on different roughness configurations. Piles with idealised triangular asperities installed in a soft rock were created, tested and imaged using X-ray CT under increased vertical loads. The X-ray CT images at different pile head displacements were qualitatively analysed and the relevant pile-rock interface load-bearing mechanisms were established. After ascertaining the relevance of the experimental shaft resistance values with the different empirical correlations, the laboratory pile-load tests were modelled using DEM to further investigate the micro-mechanic behaviour of the piles. The major conclusions are as follows: 1. The mechanisms happening at the pile-rock interface during the vertical displacement of the pile into the soft rock can be characterised into three stages based on the load resistance attributes and 3D image visualisation: Stage I-sliding; Stage II-local shearing; Stage IIIprogressive shearing. For piles with minimal interface roughness, the mechanism is dominated by sliding, whereas for the rougher piles, it is dominated by asperity failure. At the end of these stages, the shaft carries a residual load, and its magnitude depends upon the roughness of the rock socket. 2. In modelling the experimental load test of piles using DEM, the bonded particle model failed to capture the mechanical behaviour of the soft rocks due to the absence of softening micro-parameters that control the damage evolution in the post-peak region. However, the cohesive damage plasticity model can be used to simulate the rock media effectively as it can capture the macro-response of soft rocks with reasonable accuracy. 3. The effect of the symmetrical DEM model sections on the pile behaviour was investigated. It can be inferred that using a 45° section overestimates the total resistance by 30%, whereas adopting a 90° section or higher significantly reduces the boundary effect.
4. The load-displacement behaviour at the shaft is predominantly controlled by the leading edges of the pile asperity when compared to the trailing edges. The horizontal and vertical forces acting on these segments increase with the depth of the pile. The bottom leading edge segment (S6), in particular, resists a substantial portion of the shaft load through a continuous strain-hardening approach. It can be attributed to its connectivity with the base rock and the absence of an opening to generate a shear plane. 5. The expansion of the cylindrical cavity of the socket predominantly falls in the elastic zone for piles with low roughness (7.5°) and in the plastic zone for higher roughness profiles (17.5°, 22.5°). The height of the asperity determines the amount of debris expelled at the Fig. 18 Micro-crack evolution from DEM: a count for different roughness profiles, b development of different failure modes in 17.5° roughness configuration for each stage of the loadbearing mechanism pile-rock interface during loading, whose compression ultimately transforms as the residual resistance at the shaft. The piles with very low roughness will experience a compression force opposite to pile loading due to the elastic rebound of the rock at the pile-rock interface, which gains significance when the vertical pile displace-ment is more than half asperity length. Though the effect is minimal, it results in the reduction of the shaft capacity at large displacements. However, for the piles with higher roughness, this phenomenon is not observed as the asperities are sheared down into the shaft as debris. 6. The DEM results concerning the micro-crack formation and damage evolution at the interface support the mechanisms interpreted from the X-ray CT images. The asperity zones depict an upsurge in the number of yielding contacts with an increase in the asperity height. It was interpreted from the micro-crack count that shearing forms the primary basis of failure at the pile-rock interface, followed by the mixed-mode failure. Also, the piles with low roughness (7.5°), though yielding at higher loads, tend to behave like the smooth piles (0°), as observed from the micro-crack evolution at the interface.
These insightful micro-mechanical contributions at the asperity level are believed to improve the understanding of the pile-rock interactions. Additional experimental and numerical studies may be required to explore for any possible geometric scale effects. Through this type of research, the knowledge on the interface mechanisms at the grainscale is enhanced, which can subsequently aid in the development of refined design criteria with the probable inclusion of construction-induced defects (e.g. smear fabric).
Funding Open Access funding enabled and organized by CAUL and its Member Institutions. This study was funded through the Monash International Tuition Scholarship (MITS) and the Research Training Program (RTP) stipend provided by the Australian Government as a part of the PhD programme at Monash University.

Data Availability
The data supporting the findings of this study are available from the corresponding author upon reasonable request.

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