## 2 Proposed \"Hybrid BEM-p-y Curve\" Method
### Non-Linear Solution Procedure  
The solution system is defined as: [F][X] = [P] (Equation (5)). [X] is the unknowns vector composed of \(k+2\) or \(k+1\) terms for free or fixed head conditions, respectively, where \(k\) is the number of pile blocks,  
Figure 4: Mindlin solution scheme. Figure 5: Hybrid BEM-p-y method. \(p\) are the \(k\) unknown pressures acting at the pile-soil interface, \(y_{0}\) is the single pile displacement at the pile-head, \(\theta_{0}\) is the pile-head rotation. [P] is the known term vector, which has the same dimension as the vector [X]. [F] is a \((k+2)\)\(x\)\((k+2)\) or \((k+1)\)\(x\)\((k+1)\) matrix, obtained by adding the \(k\)\(x\)\(k\) pile flexibility matrix [F\({}_{\text{P}}\)] (containing the \(a_{ij}\) coefficients), the \(k\)\(x\)\(k\) soil (\"far field\") flexibility matrix [F\({}_{\text{SFF}}\)] (composed of the \(b_{ij}\) coefficients), and the \(k\)\(x\)\(k\) soil \"near field\" flexibility matrix [F\({}_{\text{SNF}}\)] (composed of the \(c_{ii}\) coefficients, representing the flexibility of each non-linear spring). The final two (or one in the fixed-head condition) rows and columns of the matrix [F] are necessary to impose the equilibrium and to complete the compatibility equations at each pile-spring-soil interface node, respectively. \[\left[ \begin{array}{cccc}a_{11}+b_{11}+c_{11}&a_{1j}+b_{1j}&a_{1k}+b_{1k}&- 1&-z_{i}\\ a_{i1}+b_{i1}&a_{ii}+b_{ii}+c_{ii}&a_{ik}+b_{ik}&-1&-z_{i}\\ a_{k1}+b_{k1}&a_{kj}+b_{kj}&a_{kk}+b_{kk}+c_{kk}&-1&-z_{k}\\ 1&1&1&0&0\\ z_{1}&z_{j}&z_{k}&0&0\end{array} \right]\left[ \begin{array}{c}p_{1}\cdot \Delta_{1}D\\ p_{i}\cdot\Delta_{i}D\\ p_{k}\cdot\Delta_{k}D\\ y_{0}\\ \theta_{0}\end{array} \right]=\left[ \begin{array}{c}0\\ 0\\ 0\\ H\\ M\end{array} \right] \tag{5}\]  
In Equation (5), \(H\) is the horizontal load applied and \(M\) is the bending moment applied. The [F\({}_{\text{SNF}}\)] is a diagonal matrix, and each coefficient represents the \"p-y\" tangent flexibility evaluated at each spring displacement \"y\" value reached at each pile-spring node at the previous load increment (or sub-increment) step. Elements of the soil \"far-field\" matrix [F\({}_{\text{SFF}}\)] always remain constant, and are used to consider the interaction between the non-linear springs. The pile flexibility matrix, [F\({}_{\text{P}}\)], only in the case of a non-linear \"moment-curvature\" relationship for the pile section, is updated at each step, using the tangent flexural rigidities of the section, according to the bending moments reached at each pile node in the previous load increment (or sub-increment). In addition, at each load step (or sub-step) a check is carried out to determine whether the ultimate soil resistance at the pile-soil interface has been reached. Once the initial flexibility matrix has been calculated, the total horizontal load is applied in the first step of the solution. At each generic load increment \(h_{k}\), an iterative process is performed where two solutions are obtained, the first using \(h_{k}\) as the load increment, the second using two load steps equal to \(h_{k}/2\). The iterative scheme is described in Figure 6, which, for the sake of simplicity refers to the explicit Euler method with step-doubling and adaptive step-size control. However, a fourth order Runge-Kutta method can also be used to obtain some improvement in the accuracy of the solution. Once these two solutions have been computed, the incremental ratio, \(\varepsilon\), is computed according to Equation (6). \[\varepsilon=\frac{\Delta u_{2}-\Delta u_{1}}{\Delta u_{1}} \tag{6}\]  
where \(\Delta u_{1}\) and \(\Delta u_{2}\) are the incremental displacement at the pile-head evaluated using one and two steps, respectively. The \(\varepsilon\) value is compared with a predefined tolerance taken as equal to 0.001 (Figure 7).