## 2 Model Formulation and Analysis
### Mathematical Modeling  
The solution proposed by [PERSON] et al. [34] is important for providing the settlement profile at the ground surface of the column-reinforced foundation. However, the limitation is that this solution relies on the assumed deformation shape function, which was constructed using simple functions. Therefore, the solution does not satisfy the stress equilibrium condition of the model. To overcome this limitation, an improved solution is proposed by using the variational principles, for which the deformation shape is obtained naturally. In addition, the equilibrium condition is also ensured by the variational principles. In this study, the soil is assumed to behave as a linearly elastic homogeneous continuum, with a constant modulus \(E_{s}\) and a constant Poisson's ratio \(\
u_{s}\), unaffected by the presence of columns, and remaining unchanged throughout the loading process. It is also assumed that the column material behaves as a linearly elastic homogeneous material with a deformation modulus \(E_{c}\) and a constant Poisson's ratio \(\
u_{c}\). Furthermore, it is assumed that the interface between soil and column remains intact and does not slip, and the displacements of the column and the soil at the interface are equal. Since the unit cell is symmetric, only half of the model is considered in the analysis, as shown in Figure 4 in which the depth is equal to \(H\) and the radius of the column is \(r_{c}\). The radius of the entire model is \(d_{e}/2\), which is determined by Equation (1) with the spacing of the columns, \(s\), based on the area equivalence assumption. When the stress, \(\sigma_{i}\), is applied to the column, and the stress, \(\sigma_{o}\), is applied to the surrounding soil, the column-reinforced soil starts to deform. As linear elastic behavior is assumed, the total potential energy \(\Pi\) of the entire system can be introduced as  
\[\Pi=\prod_{c}+\int\limits_{r_{c}}^{\frac{d_{e}}{2}}\int\limits_{0}^{\frac{d_{ e}}{2}}\int\limits_{0}^{H}\frac{1}{2}(\sigma_{e}\varepsilon_{z}+\sigma_{r} \varepsilon_{r}+\tau_{rz}\gamma_{rz})rdrdzd\theta-\int\limits_{0}^{2\pi}\int \limits_{r_{c}}^{\frac{d_{e}}{2}}\sigma_{o}W(r,0)rdrd\theta \tag{3}\]  
Figure 4: Axisymmetric model of the column-reinforced foundation. where \(\Pi_{c}\) is the total energy of column inclusion, including the internal and external energy, and \(\sigma_{z}\), \(\sigma_{r}\), \(\tau_{rz}\), \(\epsilon_{z}\), \(\epsilon_{r}\), \(\gamma_{rz}\) are the stresses and corresponding strains of a small element of the soil continuum. \(W(r,0)\) is the vertical displacement at the surface of the surrounding soil at the location of \((r,0)\) in the cylindrical coordinate system. The constitutive equations in an axisymmetric condition can be written as  
\[ \begin{pmatrix}\sigma_{\rm r}\\ \sigma_{\theta}\\ \sigma_{z}\\ \tau_{rz}\end{pmatrix} =\frac{(1-\
u)E}{(1+\
u)(1-2\
u)} \begin{bmatrix}1& \frac{\
u}{1-\
u}&\frac{\
u}{1-\
u}&0\\ \frac{\
u}{1-\
u}&1&\frac{\
u}{1-\
u}&0\\ \frac{\
u}{1-\
u}&\frac{\
u}{1-\
u}&1&0\\ 0&0&0&\frac{1-2\
u}{2(1-\
u)}\end{bmatrix}\begin{pmatrix}\epsilon_{r}\\ \epsilon_{\theta}\\ \epsilon_{z}\\ \tau_{rz}\end{pmatrix}  \tag{4}\]  
Because the peripheral displacement equals zero due to the axisymmetric condition, the corresponding strain can be described as  
\[\epsilon_{r}=\frac{\partial u(r,z)}{\partial r},\epsilon_{\theta}=\frac{u(r,z )}{r},\epsilon_{z}=\frac{\partial W(r,z)}{\partial z},\gamma_{rz}=\frac{ \partial u(r,z)}{\partial z}+\frac{\partial W(r,z)}{\partial r} \tag{5}\]  
where \(u(r,z)\) and \(W(r,z)\) represent the radial displacement and the vertical displacement, respectively. To simplify the analysis, the following assumptions are adopted [37]. 1. The radial displacement \(u\) is assumed to be equal to zero everywhere throughout the soil column model. Because the radial displacement is negligible compared to the vertical displacement under vertical load conditions. 2.