# A macro-element for circular shallow foundations on rigid inclusion-reinforced soil

Ramon Alcala-Ochoa, Zheng Li, Panagiotis Kotronis, Giulio Sciarra

► To cite this version:

Ramon Alcala-Ochoa, Zheng Li, Panagiotis Kotronis, Giulio Sciarra. A macro-element for circular shallow foundations on rigid inclusion-reinforced soil. Soil Dynamics and Earthquake Engineering, 2026, 202, pp.110037. 10.1016/j.soildyn.2025.110037. hal-05454871

HAL Id: hal-05454871

https://cnrs.hal.science/hal-05454871v1

Submitted on 12 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A macro-element for circular shallow foundations on rigid inclusion-reinforced soil

Ramon Alcala-Ochoa<sup>a</sup>, Zheng Li<sup>b,D</sup>, Panagiotis Kotronis<sup>a,D,*</sup>, Giulio Sciarra<sup>a,D</sup>

## ARTICLE INFO

### Keywords:

- Macro-element

- Rigid inclusions

- Shallow foundations

- Reinforced soils

- Ground improvement

## ABSTRACT

This study presents a plasticity-based macro-element (ME) for circular shallow foundations on rigid-inclusion (RI)-reinforced, layered soils subjected to seismic loading. The formulation includes: (i) an inclined interaction ellipse in the (h, m) plane with closed-form sizing and tilt (a, b, ψ) linked to RI descriptors (α, $H_{LTP}$, $φ_{LTP}$); (ii) elastic terms extracted from dynamic impedance functions sampled at the soil-structure interaction frequency; and (iii) a vertical hardening variable γ whose evolution included contributions from plastic horizontal displacement and rotation. A diagonal elastic matrix was adopted as a first-order approximation for shallow circular footings with small eccentricities. The model was validated against dynamic centrifuge tests on RI-reinforced ground: moment-rotation hysteresis and cyclic degradation were reproduced consistently, top-displacement histories were matched to a reasonable degree, and settlement trends were captured with a conservative response. The failure-envelope fit is limited to circular footings $α \in [2.5, 7.0]\%$, $H_{LTP} \in [0.5, 1.0]$ m, and $φ_{LTP} \in [30°, 42°]$. Within this domain, the ME offered an efficient tool for nonlinear soil-structure interaction analysis of RI-supported systems.

## 1. Introduction

Rigid inclusions (RI) are stiff elements, typically in the form of columns or piles, not connected with the foundation, introduced into soft soil to enhance its mechanical properties. Beyond the clear benefit of reduced settlement (vertical displacement), the proper integration of RI into the foundation design mitigates the risk of collapse, enhances the overall stability and improves the structural seismic response. In numerous cases, RI foundations have proven to be a more efficient and cost-effective alternative to traditional connected pile foundations [1-4].

Currently, numerical simulations of RI foundations under various loading conditions using the Finite Element Method (FEM) or the Finite Difference Method (FDM) are widely adopted within research and engineering communities. For instance, Saadatinezhad et al. [5] conducted a parametric study via FEM simulations on various lengths, diameters, spacings, and numbers of RI to investigate their effects on the load transfer platform and the corresponding load transfer mechanisms. Similarly, Mánica Malcom et al. [6] performed numerical parametric studies using FDM to analyse ground responses due to a group of RI in soft Mexico City clay. Numerous similar numerical studies can be

found in the literature for RI foundations under both static and dynamic loadings [7-11].

However, FEM or FDM simulations consume significant computational time and resources. Under complex loading conditions, evaluating soil-structure interaction (SSI) phenomena through FEM or FDM simulations remains challenging. As an alternative solution, macro-element (ME) modelling has shown promising potential in addressing this issue. The concept of ME was originally proposed by Nova and Montrasio [12] to consider soil-shallow foundation interactions through a 3D upscaling constitutive law in the generalized force-displacement space. In this paper, the term generalized failure envelope denotes an upper-scaled system-level limit surface that bounds admissible combinations for the foundation-soil system in the space of generalized forces (v, h, m); it does not correspond to a conventional soil yield surface.

In their pioneering work, Nova and Montrasio [12] lumped the nonlinear SSI phenomena of a shallow foundation and the surrounding soil in a single point at the base of the structure. Consequently, the computation time for SSI problems is reduced by thousands of times compared to conventional FEM or FDM simulations. ME can be applied

<footnote><sup>a</sup> Nantes Université, Ecole Centrale Nantes, Institut de Recherche en Génie Civil et Mécanique (GEM), CNRS, UMR-6183, 1 Rue de la Noë, Nantes, 44000, France</footnote>

<footnote><sup>b</sup> Laboratoire Centrifugeuses Géotechniques CG, Département GERS, Université Gustave Eiffel Campus de Nantes (ex IFSTTAR), Allée des Ponts et Chaussées, 44340 Bouguenais, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: ramon.alcala-ochoa@outlook.com (R. Alcala-Ochoa), zheng.li@univ-eiffel.fr (Z. Li), panagiotis.kotronis@ec-nantes.fr (P. Kotronis), giulio-sciarra@ec-nantes.fr (G. Sciarra).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. The macro-element concept for a RI foundation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Generalized forces and generalized displacements at the footing of a RI foundation, two-dimensional problem.</figcaption></figure>

for monotonic, cyclic and dynamic loadings. Due to its high computational efficiency, ME approach provides accurate and rapid assessments of SSI problem response, enabling fast design and decision-making for engineers.

In the past 30 years, the ME approach has been developed and applied with great success in simulating SSI problems of shallow foundations [13–23], offshore plate-anchors and spud-can foundations [24, 25], caisson foundations [26], pipelines [27], piled foundations [11, 28–33], but also structural elements (e.g., reinforced concrete beam and column sections [34]).

This study aims to develop a novel ME to simulate the nonlinear behaviour of RI foundations subjected to static and seismic loadings. The ME presented in this paper is based on the mathematical theory of plasticity [35,36], the ME concept introduced in geomechanics by Nova and Montrasio [12], and the existing macro-elements for shallow foundations developed by Cremer et al. [14,15], Grange [17] Grange et al. [18,19,20,37] and Salciarini et al. [23].

The mathematical formulation of the elasto-plastic ME is presented in Section 2. This is followed by an extension to dynamic (seismic loading) problems in Section 3, a synopsis of all the formulas in Section 4, and the methods for calibrating the ME parameters in Section 5. Finally, the performance of the new ME is validated using dynamic centrifuge test data [38,39] in Section 6.

## 2. Mathematical formulation

According to the ME concept, the nonlinear behaviour of a shallow circular RI foundation, the load transfer platform (LTP), and the surrounding soil is concentrated at a single point at the base of the structure, as shown in Fig. 1.

For a two-dimensional problem the generalized force **F** and the generalized displacement **u** vectors are written as follows:

$$ \underline{\underline{\mathbf{F}}} = \begin{bmatrix} \underline{v} \\ \underline{h} \\ \underline{m} \end{bmatrix} = \frac{1}{V_0} \begin{bmatrix} \underline{V} \\ \underline{H} \\ \underline{M/B} \end{bmatrix} \qquad \underline{\underline{\mathbf{u}}} = \begin{bmatrix} \underline{w} \\ \underline{u} \\ \underline{\theta} \end{bmatrix} = \frac{1}{B} \begin{bmatrix} \underline{u}_z \\ \underline{u}_x \\ \underline{B\theta}_y \end{bmatrix} \qquad (1) $$

with **V** the vertical force, **H** the horizontal force, **M** the moment, **u**_z the vertical displacement, **u**_x the horizontal displacement and **θ**_y the rotation respectively, see Fig. 2. In Eq. (1), $V_0$ is the vertical bearing capacity, **B** the diameter of the circular foundation, **v**, **h**, **m** the dimensionless forces and **w**, **u**, **θ** the dimensionless displacements at the centre of the footing. Note that with this type of normalization,

**F** and **u** are conjugate; the work of the external forces applied to the foundation can be expressed using these variables [12].

Following the mathematical theory of plasticity [35,36], the generalized displacement vector **u** can be decomposed in an elastic **u**^el and a plastic part **u**^pl as follows:

$$ \underline{\underline{\mathbf{u}}} = \underline{\underline{\mathbf{u}}}^{\text{el}} + \underline{\underline{\mathbf{u}}}^{\text{pl}} \qquad (2) $$

### 2.1. Elastic behaviour

The relationship between the load **F** and displacement **u**^el in the elastic range can be written as follows:

$$ \underline{\underline{\mathbf{F}}} = \underline{\underline{\mathbf{K}}}^{\text{el}} \underline{\underline{\mathbf{u}}}^{\text{el}} = \begin{bmatrix} \underline{v} \\ \underline{h} \\ \underline{m} \end{bmatrix} = \begin{bmatrix} k_1 & 0 & 0 \\ 0 & k_2 & 0 \\ 0 & 0 & k_3 \end{bmatrix} \begin{bmatrix} \underline{w}^{\text{el}} \\ \underline{u}^{\text{el}} \\ \underline{\theta}^{\text{el}} \end{bmatrix} \qquad (3) $$

where **K**^el is the elastic stiffness matrix and the diagonal terms $k_1$, $k_2$, $k_3$ the stiffness associated with the vertical, horizontal and rotational motion, respectively. For a shallow, surface-founded circular footing subjected to combined (V, H, M), the translation-rotation coupling terms in **K**^el are typically second order compared with the diagonals [14,17,19,21,22].

On the present configuration, a diagonal **K**^el was adopted as in (3), which simplified calibration without degrading the predicted global response.

The diagonal assumption may be insufficient when: (i) the foundation is embedded or skirted, (ii) the base is markedly non-circular or large load eccentricities prevail, (iii) layering/anisotropy or sloping ground induces static coupling, or (iv) the response is governed by high-frequency rocking close to the SSI mode. In such cases, the off-diagonal terms can be obtained from small-amplitude FE impedance analyses and inserted directly into **K**^el.

The rate form of Eq. (3) is written as (where the symbol [·] represents the time derivative):

$$ \begin{bmatrix} \dot{\underline{v}} \\ \underline{h} \\ \dot{\underline{m}} \end{bmatrix} = \begin{bmatrix} k_1 & 0 & 0 \\ 0 & k_2 & 0 \\ 0 & 0 & k_3 \end{bmatrix} \begin{bmatrix} \dot{\underline{u}}^{\text{el}} \\ \dot{\underline{u}}^{\text{el}} \\ \dot{\underline{\theta}}^{\text{el}} \end{bmatrix} \qquad (4) $$

where the dimensionless forms of the stiffness parameters are given as:

$$ \left\{ \begin{aligned} k_1 &= \frac{k_z^{\text{el}} B}{V_0} \\ k_2 &= \frac{k_x^{\text{el}} B}{V_0} \\ k_3 &= \frac{k_{\theta}^{\text{el}}}{B V_0} \end{aligned} \right. \qquad (5) $$

The horizontal elastic stiffness $k_x^{\text{el}}$ is here approximated using the classical expression of Gazetas [40] for a circular foundation on a homogeneous elastic half-space:

$$ k_x^{\text{el}} = \frac{4G_0 B}{2 - \nu} \qquad (6) $$
where $G_0$ is the equivalent shear modulus of a homogeneous soil and $v$ is the Poisson ratio.

This simplification follows common practice in ME formulations for shallow foundations [13,17] and provides a consistent first-order estimate of the foundation-soil stiffness. In RI systems, however, the ground is stratified (load transfer platform over soft soil over stiff bearing layer), and such a configuration cannot be rigorously represented by Eq. (6). The validity and limitation of this approximation is therefore examined using the impedance results in Fig. 6. The figure shows that the horizontal dynamic stiffness obtained from the FEM model for the RI-foundation is nearly identical to that of an equivalent homogeneous soil without inclusions, suggesting that, for this configuration-where the soft soil layer governs the overall response-the inclusions have little effect on horizontal rigidity.

In contrast, the vertical and rotational stiffnesses (Fig. 6(a) and (c)) display marked differences between the homogeneous shallow and RI foundation cases even in an homogeneous soil, reflecting the increased composite stiffness induced by the inclusions. This confirms that Eq. (6) hypothesis is acceptable only for the horizontal direction in this preliminary implementation, whereas the other stiffness components must be determined from case-specific impedance analyses. The extension of this approach to layered or RI-reinforced soils could be addressed in future work through multilayer impedance formulations or dedicated 3D numerical parametric studies.

Currently, no closed-form solutions exist for evaluating the vertical and rotational elastic stiffnesses ($k_z$, $k_θ$) of RI foundations. A practical approach is to identify these parameters from linear elastic FEM simulations, by applying small imposed displacements or rotations at the foundation and measuring the resulting reaction forces. Alternatively, they can be estimated from dynamic impedance calculations, as detailed in Section 3. For preliminary assessments on RI-reinforced or layered soils, the horizontal stiffness $k_x^{el}$ can be taken from the Gazetas expression in (6), whereas the vertical and rocking stiffnesses $k_z^{el}$ and $k_θ^{el}$ must be identified from small-amplitude FE impedance analyses (or from project-specific charts derived therefrom).

## 2.2. Plastic behaviour

Following the plasticity theory [35,36], the plastic components are expressed as:

$$ \begin{bmatrix} \ddot{u}^{pl} \\ \ddot{u}^{pl} \\ \dot{\theta}^{pl} \end{bmatrix} = \lambda \begin{bmatrix} \frac{\partial g}{\partial \nu} \\ \frac{\partial g}{\partial h} \\ \frac{\partial g}{\partial m} \end{bmatrix} \quad (7) $$

where $\lambda$ represents the plastic multiplier, determining the magnitude of the increment of plastic displacements, while the derivatives of the plastic potential $g$ define its direction.

The key components required to formulate the plastic behaviour of the ME include a failure/ultimate bearing capacity envelope, the loading surface and its evolution (to ensure a smooth transition after yielding to the failure envelope) and the plastic flow rule. These elements are detailed in the following sections.

### 2.2.1. Failure envelope

The failure envelope discussed herein is the one introduced by Alcala-Ochoa et al. [41,42] and it is briefly described below. In the context of ME modelling, this generalized failure envelope does not represent a conventional soil failure surface. Instead, it corresponds to a system-level, force-based limit surface that defines the boundary of admissible generalized forces for the foundation-soil system.

If $\alpha$ is the coverage area, $H_{LTP}$ the LTP thickness and $\varphi_{LTP}$ the LTP friction angle, the 3D failure envelope of a RI foundation takes the form:

$$ f_{\infty} := f(v, m, h) = \frac{h^2}{h_0^2} + \frac{m^2}{m_0^2} - 2a_0 \frac{hm}{h_0 m_0} - [4v(1 - v)]^2 = 0 \quad (8) $$

<table><caption>**Table 1**
Fitting parameters of $h_0$, $m_0$ and $a_0$.</caption><thead><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>a<sub>4</sub></th></tr></thead><tbody><tr><th>h<sub>0</sub></th><td>0.1717</td><td>5.083</td><td>0.5183</td><td>-5.057</td></tr><tr><th>m<sub>0</sub></th><td>0.1073</td><td>10.45</td><td>0.4129</td><td>-4.284</td></tr><tr><th>a<sub>0</sub></th><td>-0.3952</td><td>-31.06</td><td>-0.7157</td><td>6.247</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Normalized failure envelope in $v - m - h$ space. Proposed analytical formula (yellow, Eq. (8)) compared with numerical sweep test results (grey dots [41,42]).</figcaption></figure>

where the fitting parameters of $h_0$, $m_0$ and $a_0$ in Eq. (8) are estimated by:

$$ c = a_1 \left[ 1 + a_2 \left( \frac{\alpha}{H_{LTP}} \right)^2 + a_3 (\varphi_{LTP})^2 + a_4 \left( \frac{\alpha}{H_{LTP}} \varphi_{LTP} \right) \right] \quad (9) $$

Here $c$ represents the parameters $h_0$, $m_0$ and $a_0$, and $a_1$, $a_2$, $a_3$, $a_4$ are summarized in Table 1:

The normalized failure envelope in the $v - m - h$ space is shown in Fig. 3.

The main hypotheses and validity domain of Eq. (8), Eq. (9) and Table 1 are summarized as follows:

* A shallow circular rigid foundation.

* Linear elastic behaviour for the shallow foundation and the RI.

* A three-layer stratigraphy (LTP + soft soil + stiff soil) with specific material properties, as described in [41].

* A coverage area $\alpha$ ranging from 2.5% to 7.0%.

* An LTP thickness $H_{LTP}$ ranging from 0.5 m to 1.0 m.

* An LTP friction angle $\varphi_{LTP}$ ranging from 30° to 42°.

* The maximum size ellipse is calibrated at $v = 0.5$, as detailed by Alcala-Ochoa [41].

* The failure surface is determined for the maximum bearing capacity, $V_0$, defined using the tangent method. It is supposed to be consistent for the other force levels, as discussed by Alcala-Ochoa [41].

Eq. (8) can be rewritten as:

$$ f_{\infty} = \left( \frac{h}{4h_0 v(1 - v)} \right)^2 + \left( \frac{m}{4m_0 v(1 - v)} \right)^2 - 2a_0 \frac{hm}{16h_0 m_0 v^2 (1 - v)^2} - 1 = 0 \quad (10) $$

This adheres to the general form

$$ f_{\infty} = \left( \frac{h}{av^c(1 - v)^d} \right)^2 + \left( \frac{m}{bv^e(1 - v)^f} \right)^2 + \psi \frac{hm}{abv^{(c+e)}(1 - v)^{(d+f)}} - 1 = 0 \quad (11) $$

where $a$ and $b$ define the size of the ellipse and $\psi$ its inclination in the $m - h$ space; $c$, $d$ and $e$, $f$ control the parabolic shape in the $v - h$ and $v - m$ plane, respectively. To ensure that Eq. (10) matches Eq. (11), the parameters should be set as $c = d = e = f = 1$.

Eq. (11) is therefore similar to the “rugby ball shape” failure envelope, originally attributed to Nova and Montrasio [12] and employed by researchers in the literature [13,17,43–45]. The novelty here is that Eq. (11), designed for RI foundations, accounts for the inclination of the ellipse in the $m - h$ space through the parameter $\psi$.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Eq. (13) Exponential relation and FE results for different normalized vertical forces for a circular RI foundation.</figcaption></figure>

### 2.2.2. Loading surface

Similar to the macro-elements for shallow foundations proposed by Cremer [13], Cremer et al. [14,15], Grange [17] and Grange et al. [18,19,20], the loading surface displays a shape akin to the failure envelope of Eq. (11). Its progression is characterized by a series of hardening variables as follows:

$$f_c(\underline{\underline{F}}, \underline{\underline{\tau}}, \rho, \gamma) := \left( \frac{h}{\rho a v^c (\gamma - \upsilon)^d} - \frac{\delta}{\rho} \right)^2 + \left( \frac{m}{\rho b v^e (\gamma - \upsilon)^f} - \frac{\beta}{\rho} \right)^2 + \psi \left( \frac{h}{\rho a v^c (\gamma - \upsilon)^d} - \frac{\delta}{\rho} \right) \left( \frac{m}{\rho b v^e (\gamma - \upsilon)^f} - \frac{\beta}{\rho} \right) - 1 = 0 \quad (12)$$

where $\underline{\underline{\tau}} = (\delta, \beta)$ is a vector containing the two kinematic variables $\delta$ and $\beta$, $\rho$ is an isotropic hardening variable and $\gamma$ is a variable that describes the expansion of the loading surface along the $\upsilon$ direction.

#### Kinematic hardening variables $\beta$ and $\delta$

The kinematic hardening variables $\beta$ and $\delta$ dictate the evolution of the ellipse's centre in the $m-h$ plane. Given the similarity in their equations, only the evolution of $\delta$ is given hereafter.

As it is noted by Grange [17], both experimental [45] and numerical results [13] for circular shallow foundations on various soil types under monotonic static loadings indicate that the force-displacement curves include an exponential term and are independent of the vertical force. Consequently, the evolution of the internal variables can be derived using first-order differential equations. For example, the relationship between the horizontal force $h$ and the horizontal displacement $u$ can be expressed as:

$$h = h_{\infty}^* \left( 1 - \exp \left( - \frac{k_2}{h_{\infty}^*} u \right) \right) \quad (13)$$

where $h_{\infty}^* (0 < h_{\infty}^* \le 1)$ represents the curve limit of $h$ when $u$ tends to infinity.

As shown by three-dimensional FEM calculations in the numerical study by Alcala-Ochoa [41], circular shallow foundations on soil exhibit similar behaviour in the horizontal direction, regardless of the presence of RI. Thus, Eq. (13) is also applicable to RI foundations. For example, Fig. 4 presents the representation of Eq. (13) using green bold lines, alongside FE numerical results for various normalized vertical forces for a circular RI foundation.

The associated 1st order differential equation has the following form:

$$\dot{h} = k_2 \left( 1 - \frac{h}{h_{\infty}^*} \right) \dot{u} \quad (14)$$

Considering Eq. (2) and that $\underline{\underline{F}} = \underline{\underline{K}}^{\text{el}}\underline{\underline{u}}^{\text{el}}$ we get:

$$\dot{h} = k_2 \left( \frac{h_{\infty}^{*}}{h} - 1 \right) \dot{u}^{\text{pl}} \quad (15)$$

The evolution of the kinematic hardening variable $\delta$ can be therefore defined as [17]:

$$\dot{\delta} = \frac{k_2}{2av^c (\gamma - \upsilon)^d} \left( \frac{h_{\infty}^*}{\delta} - 1 \right) \dot{u}^{\text{pl}} \quad (16)$$

For cyclic (quasi-static) or seismic (dynamic) loadings, the behaviour of the foundation is dependent on the loading direction. For instance, when subjected to severe loading in the right direction, the distribution of soil plastic strains varies beneath the right and left parts of the foundation. Upon load reversal, the asymmetric distribution of stresses affects the foundation's behaviour. To incorporate the influence of the loading direction, the evolution of the kinematic hardening variables can be expressed as [17]:

$$\left\{ \begin{array}{l} \dot{\delta}^{\oplus} = \frac{k_2}{2av^c (\gamma - \upsilon)^d} \left( \frac{h_{\infty}^{*}}{\delta^{\oplus}} - 1 \right) \quad |\dot{u}^{\text{pl}}| \\ \\ \dot{\delta}^{\ominus} = \frac{k_2}{2av^c (\gamma - \upsilon)^d} \left( -\frac{h_{\infty}^{*}}{\delta^{\ominus}} - 1 \right) \quad |\dot{u}^{\text{pl}}| \end{array} \right. \quad (17)$$

where at each step, $\delta$ is computed as:

$$\delta = \delta^{\oplus} + \delta^{\ominus} \quad (18)$$

$h_{\infty}^*$ is calculated as:

$$\left\{ \begin{array}{ll} h^{*\infty} = \delta_{\text{lim}} - \delta^{\ominus} & \text{if } \delta^{\ominus} > 0 \\ h^{*\infty} = -(\delta_{\text{lim}} - \delta^{\ominus}) & \text{if } \delta^{\ominus} < 0 \end{array} \right. \quad (19)$$

where $\delta_{\text{lim}}$ is determined using the tangency rule [13,46], utilized to ensure that the loading surface does not intersect the failure criterion.

#### Isotropic hardening variable $\rho$

A relationship can be introduced between the isotropic $\rho$ and the kinematic hardening variables $\underline{\underline{\tau}}$ [13,17], (where the symbol $|||$ represents the norm and $d./dt$ the symbol for the derivative with respect to time)

$$\dot{\rho} = \left| \frac{d \| \underline{\underline{\tau}} \|}{dt} \right| = \left| \frac{\delta \dot{\delta} + \beta \dot{\rho}}{\sqrt{\delta^2 + \rho^2}} \right| \quad (20)$$

that translates the fact that when a plastic state is reached during a new cycle, plastic behaviour is recovered at the same state and with the same slope.

#### Variable $\gamma$

The foundation is initially submitted only to the structure's self-weight ($h = m = 0$) and therefore $\gamma = \upsilon$. Then, the following empirical relationship provided by Nova and Montrasio [12] is adopted:

$$\gamma = 1 - e^{(-k_1 \upsilon)} \quad (21)$$

Eq. (21) can be written as a 1st order differential equation as:

$$\dot{\gamma} = k_1 (1 - \gamma) \dot{\upsilon} \quad (22)$$

or

$$\dot{\gamma} = k_1 \frac{1 - \gamma}{\gamma} \dot{u}^{\text{pl}} \quad (23)$$

Under the influence of overturning moments and horizontal loads, additional plastic displacements – such as horizontal displacements and rotations – can lead to the expansion of the loading surface in the $\upsilon$ direction. Consequently, the evolution of $\gamma$ can be articulated as follows [17]:

$$\dot{\gamma} = (c_1 k_1 \dot{u}^{\text{pl}} + c_2 k_2 |\dot{u}^{\text{pl}}| + c_3 k_3 |\dot{\theta}^{\text{pl}}|) \frac{1 - \gamma}{\gamma} \quad (24)$$

where $c_1, c_2, c_3$ to be determined, representing the contribution of each plastic displacement component.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Plastic potential g, the loading surface f_c, and the associated normal vectors for a specific loading point F in the v-h plane.</figcaption></figure>

Source: Modified from Cremer [13]

### 2.2.3. Flow rule

Consistent with the macro-elements introduced by Cremer [13] and Grange [17], an associate flow rule is used for the h-m plane, while a non-associate flow rule is utilized for the h-v and m-v planes. The plastic potential, as outlined by Grange [17], is given by:

$$g(\underline{\underline{F}}, \underline{\underline{\tau}}, \rho, v_c) := \begin{cases} \frac{h}{\rho \frac{av_c^c}{(\kappa\gamma + v_c)^d} (\gamma - v)^d (\kappa\gamma + v)^d} - \frac{\delta}{\rho} \\ + \left( \frac{m}{\rho \frac{bv_c^e}{(\xi\gamma + v_c)^f} (\gamma - v)^f (\xi\gamma + v)^f} - \frac{\beta}{\rho} \right)^2 \end{cases} - 1 = 0 \quad (25)$$

where $v_c$ is the current vertical force, representing the first coordinate of the loading point F, see Fig. 5. In the region of low v values, the outward normal can be negative regardless of the direction of the loading increment. Consequently, even as v increases, w decreases. This will cause the uplift at a low vertical loading level, which is not the case in reality. Opting for a plastic potential as the one presented in Eq. (25) allows increasing w^pl values within this zone, see Fig. 5.

The parameters $\kappa$ and $\xi$ are introduced to adjust the shape of the plastic potential, ensuring it features a horizontal tangent at a user-specified vertical reduced force v. Generally, the horizontal tangent in the m-v plane is located at $v = (\frac{1-\xi}{2})$, while in the v-h plane, it is positioned at $v = (\frac{1-\kappa}{2})$. As illustrated in Fig. 5, the representation of g in the v-h plane indicates that $\kappa$ and $\xi$ are set to 1, resulting in a consistent horizontal tangent at v = 0. Similar representations are applicable to the v-m plane.

## 3. Extension to dynamics

In earthquake engineering, impedance functions are utilized to account for dynamic SSI [47]. These functions are typically expressed in the following form:

$$Z(\omega) = K(\omega) + iC(\omega) \quad (26)$$

For time-domain use, constant parameters were extracted at the soil-structure interaction frequency $\omega_{SSI}$:

$$K_{ii} := \Re\{Z_{ii}(\omega_{SSI})\}, \quad C_{ii} := \frac{\Im\{Z_{ii}(\omega_{SSI})\}}{\omega_{SSI}}. \quad (27)$$

The real component, $K(\omega)$, is known as the dynamic stiffness and reflects the stiffness and inertia of the supporting soil. The imaginary

component, $C(\omega)$, represents the damping effect, primarily due to radiation within the soil. The frequency dependence in impedance functions arises from the influence of frequency on inertia, as soil properties are largely invariant with respect to frequency variations [48].

This approach is widely adopted in contemporary design methodologies, especially within substructural methods. In practice, implementations range from simple spring-dashpot fits calibrated at a representative frequency [11,49] to the classical “monkey-tail” rational fits-compact, lumped descriptions of the soil support [47,50,51]. More generally, these lumped-parameter models (LPMs) preserve the essential frequency dependence of stiffness and radiation damping over the SSI-relevant band while remaining easy to embed in time-domain analyses. In our framework, such representations can be used as radiation add-ons in parallel with the elasto-plastic ME: the imaginary component of the matched impedance represents radiation damping (implemented either as a single dashpot at a representative frequency or via a compact LMP), while the elasto-plastic ME captures the nonlinear, path-dependent soil-foundation behaviour.

To extend the proposed ME for dynamic analysis, it is proposed to calculate the impedance functions for various typical RI foundation configurations and to develop analytical equations for practical applications, like those presented by Gazetas [40,48] for shallow foundations. The real components of the impedance functions are used to calibrate the stiffness terms in Eq. (5), while the imaginary components serve to calibrate the dashpots, which can be directly incorporated into the ME as suggested by Grange et al. [20]. For use in time-domain analyses, the dynamic stiffnesses $K_{ii}(\omega)$ and radiation coefficients $C_{ii}(\omega)$ should be sampled at the soil-structure interaction frequency $\omega_{SSI}$ (first coupled mode or dominant input).

As an illustrative example, Fig. 6 presents the real components (dynamic stiffnesses) of the impedance functions for a circular shallow foundation with a diameter $B = 6$ m resting on a homogeneous soil. The impedance functions were computed through three-dimensional dynamic FEM simulations in the frequency domain. For each direction of interest (vertical, horizontal, and rocking), small-amplitude harmonic excitations were applied at the base of a rigid circular foundation, while the reaction forces and resulting displacements or rotations at the foundation interface were recorded. The complex impedance was then obtained as the ratio between the Fourier transforms of force and displacement:

$$K(\omega) = \frac{\hat{F}(\omega)}{\hat{u}(\omega)}, \quad (28)$$

from which the real part provides the dynamic stiffness and the imaginary part the radiation damping coefficient.

The numerical model was implemented in Cast3M [52] as a three-dimensional axisymmetric-by-sectors domain discretized with eight-node bricks (CUB8), refined towards the foundation-soil interface. The soil domain radius and depth were set to 12B and 5B to limit boundary reflections over the frequency band of interest. A wavelength criterion was enforced, $\Delta\ell \le \frac{v_s}{10f_{max}}$, with $v_s = 200$ m/s and $f_{max} = 15$ Hz, with a graded refinement in the near field. The rigid foundation had $B = 6$ m and thickness 0.5 m. RIs, when present, were linear-elastic cylinders ($E = 20$ GPa, $d = 0.6$ m, spacing $s = 2$ m) embedded in an otherwise homogeneous soil. Two cases were analysed: homogeneous soil with and without inclusions. This $v_s$ value was estimated using the empirical [53] correlation between shear modulus and mean effective stress, applied to the average stress level and relative density corresponding to the LTP-soil system in the centrifuge model described by Nohra [54]. It should be emphasized the this value represents an average approximate equivalent velocity for the composite reinforced ground rather than a directly measured value. No material damping was introduced; only radiation damping was captured through the boundaries.

Viscous absorbing boundaries were applied on the lateral and bottom faces using impedance-type dashpots consistent with the local soil properties; symmetry constraints were imposed where appropriate to
<figure><img src="image_9.png" /><figcaption>Fig. 6. The real components of the impedance functions for a circular shallow foundation with a diameter $B = 6$ m on a homogeneous soil ($V_s = 200$ m/s and $v = 0.3$) with (continuous lines) and without (non-continuous lines) RI: (a) vertical and (b) horizontal and (c) rotational stiffness. The FEM results are compared with the analytical formulas from [40].</figcaption></figure>

reduce computational cost without affecting wave propagation. Small-amplitude harmonic excitations were imposed at the footing in vertical ($F_z$), horizontal ($F_x$) and rocking ($M_y$) directions. For each direction, interface reactions and displacements ($u_z, u_x, \theta_y$) were recorded and the impedance was obtained. Sampling at $\Delta t = 1/(8 f_{max})$ ensured resolution up to 15 Hz. Frequency response functions were post-processed, into real and imaginary parts, to give dynamic stiffness and radiation damping.

The horizontal dynamic stiffness (Fig. 6b) was nearly unchanged by inclusions, indicating limited influence on $k_x$ for this homogeneous soft-soil-governed configuration and supporting the use of the Gazetas expression as a first-order estimate. In contrast, vertical and rocking components (Fig. 6a,c) increased markedly with inclusions, confirming that the inclusions significantly stiffen the system in these directions. Because the benchmark used a simplified homogeneous profile to enable direct comparison with closed-form half-space solutions, the curves in Fig. 6 should be interpreted as a theoretical reference; deviations with analytical results reflect stem from both the layered, inclusion-reinforced nature of the real soil and the modelling simplifications (e.g., neglect of inclusion-soil heterogeneity). A more rigorous impedance formulation would explicitly account for the layered stiffness distribution introduced by the inclusion network, which is left for future work.

It should be noted that the example in Fig. 6 does not correspond to the actual centrifuge configuration, as no load transfer platform (LTP) or underlying stiff bearing layer is included. Nevertheless, it provides a

useful benchmark to identify the qualitative effect of rigid inclusions. As observed, the presence of inclusions produces a substantial increase in the vertical and rotational stiffnesses, whereas their influence on the horizontal stiffness is minor for soft-soil conditions where the global response is governed by the deformability of the surrounding soil. Accordingly, while Eq. (6) remains adequate for preliminary estimates of $k_x^e$, the other stiffness components should be derived from dedicated impedance analyses that explicitly account for soil layering and inclusion geometry.

## 4. Synopsis

A synopsis of the proposed ME formulation follows:

* failure envelope:

$$f_{\infty} = \left( \frac{h}{av^c(1-v)^d} \right)^2 + \left( \frac{m}{bv^e(1-v)^f} \right)^2 + \psi \frac{hm}{abv^{(c+e)}(1-v)^{(d+f)}} - 1 = 0$$

* loading surface:

$$f_c(\mathbf{F}, \tau, \rho, \gamma) = \left( \frac{h}{\rho a v^c (\gamma - v)^d} - \frac{\delta}{\rho} \right)^2 + \left( \frac{m}{\rho b v^e (\gamma - v)^f} - \frac{\beta}{\rho} \right)^2 + \psi \left( \frac{h}{\rho a v^c (\gamma - v)^d} - \frac{\delta}{\rho} \right) \left( \frac{m}{\rho b v^e (\gamma - v)^f} - \frac{\beta}{\rho} \right) - 1 = 0$$
<table><caption>**Table 2**
Summary of the macro-element (ME) parameters (20).</caption><thead><tr><th>Parameter</th><th>Description</th><th>Group</th></tr></thead><tbody><tr><td><i>B</i></td><td>Diameter</td><td>Geometry</td></tr><tr><td><i>V</i><sub>0</sub></td><td>Vertical bearing capacity</td><td rowspan="4">Elastic behaviour</td></tr><tr><td><i>k</i><sub>e</sub><sup>el</sup></td><td>Vertical stiffness</td></tr><tr><td><i>k</i><sub>s</sub><sup>el</sup></td><td>Horizontal stiffness</td></tr><tr><td><i>k</i><sub>θ</sub><sup>el</sup></td><td>Rotational stiffness</td></tr><tr><td><i>C</i><sub>z</sub></td><td>Vertical dashpot</td><td>Dynamics</td></tr><tr><td><i>C</i><sub>x</sub></td><td>Horizontal dashpot</td><td>Dynamics</td></tr><tr><td><i>C</i><sub>θ</sub></td><td>Rotational dashpot</td><td>Dynamics</td></tr><tr><td><i>a</i>, <i>b</i></td><td>Shape in the <i>h</i>–<i>m</i> plane</td><td rowspan="3">Failure envelope</td></tr><tr><td><i>c</i>, <i>d</i> and <i>e</i>, <i>f</i></td><td>Shape in the <i>v</i>–<i>h</i> and <i>v</i>–<i>m</i> plane</td></tr><tr><td><i>ψ</i></td><td>Inclination in the <i>h</i>–<i>m</i> plane</td></tr><tr><td><i>c</i><sub>1</sub>, <i>c</i><sub>2</sub>, <i>c</i><sub>3</sub></td><td>Constant related to isotropic vertical variable <i>γ</i></td><td>Hardening laws</td></tr><tr><td><i>κ</i> and <i>ξ</i></td><td>Constants controlling the flow rule</td><td>Flow rule</td></tr></tbody></table>

* hardening variables:

- kinematic variables

related to the horizontal force:

$$ \begin{gathered} \delta = \delta^{\oplus} + \delta^{\ominus} \\ \dot{\delta}^{\oplus} = \frac{k_2}{2av^c(1-\nu)d} \left( \frac{h^*}{\delta^{\oplus}} - 1 \right) |\dot{u}^{pl}| \text{ and } \dot{\delta}^{\ominus} = \frac{k_2}{2av^c(1-\nu)d} \left( -\frac{h^*}{\delta^{\ominus}} - 1 \right) |\dot{u}^{pl}| \end{gathered} $$

related to the moment:

$$ \begin{gathered} \dot{\beta} = \beta^{\oplus} + \beta^{\ominus} \\ \dot{\beta}^{\oplus} = \frac{k_3}{2bt^c(1-\nu)f} \left( \frac{m^*}{\beta^{\oplus}} - 1 \right) |\dot{\theta}^{pl}| \text{ and } \dot{\beta}^{\ominus} = \frac{k_3}{2bt^c(1-\nu)f} \left( -\frac{m^*}{\beta^{\ominus}} - 1 \right) |\dot{\theta}^{pl}| \end{gathered} $$

- isotropic variable $\rho$

$$ \dot{\rho} = \left| \frac{d \| \boldsymbol{\tau} \|}{dt} \right| = \left| \frac{\delta \dot{\delta} + \beta \dot{\rho}}{\sqrt{\delta^2 + \beta^2}} \right| $$

- variable $\gamma$ related to the vertical force in the self-weight initialization:

$$ \gamma = 1 - e^{(-k_1 w)} \text{ and} $$

related to the expansion of the loading surface in $v$ direction:

$$ \dot{\gamma} = (c_1 k_1 \dot{u}^{pl} + c_2 k_2 |\dot{u}^{pl}| + c_3 k_3 |\dot{\theta}^{pl}|) \frac{1-\gamma}{\gamma} $$

* plastic potential :

$$ g(\mathbf{F}, \boldsymbol{\tau}, \rho, v_c) := \left( \frac{h}{\rho \frac{av^c}{(\kappa+v_c)d} v^c (1-\nu)d (\kappa+v_c)d} - \frac{\delta}{\rho} \right)^2 + \left( \frac{m}{\rho \frac{b v^c}{(\xi+v_c)f} v^c (1-\nu)f (\xi+v_c)f} - \frac{\beta}{\rho} \right)^2 - 1 = 0 $$

The 20 parameters of the ME, along with their respective roles, are summarized in Table 2.

## 5. Calibration

The calibration of the ME parameters is detailed below, with both simple and more complex calibration methods proposed.

(i) maximum vertical bearing capacity $V_0$: It is estimated following the semi-empirical ASIRI [55] design recommendations [55] methodology for RI-reinforced soils. In this approach, the LTP redistributes the applied load between the inclusions and the surrounding soil, and the ultimate bearing capacity is expressed as:

$$ V_0 = A_R q_R + (1 - A_R) q_S \qquad (29) $$

where $A_R$ is the inclusion area ratio, $q_R$ is the ultimate bearing stress mobilized at the inclusion tip, and $q_S$ is the bearing stress of the surrounding soft soil beneath the LTP. The value of $q_R$ is computed from classical end-bearing formulas using the soil strength parameters of the bearing layer, while $q_S$ is derived

from the bearing capacity of the soft matrix soil under the LTP. For the ASIRIplus-SDS C09 configuration ($\alpha = 7\%$, $H_{LTP} = 1$ m, $D_{RI} = 0.25$ m), this yields a prototype-scale $V_0 \approx 4.6 \times 10^6$ N, in close agreement with the value obtained from the 3D FEM tangent-method analysis. A more detailed derivation of this expression is given in [42].

(ii) elastic stiffnesses: If the stratigraphy is homogeneous, $k_x^{el}$ can be approximated with the Gazetas expression for circular foundations (Eq. (6)). In the application case, $k_x^{el}$, $k_z^{el}$, and $k_\theta^{el}$ are identified from FEM impedance analyses by reading $\Re\{Z_{ii}(\omega_{SSI})\}$. Sampling at $\omega_{SSI}$ ensured consistency with the operative SSI stiffness used in the analyses.

(iii) dashpots: Currently, there are no closed-form solutions available for RI foundations. FEM can be employed to calculate the impedances. The selection of damping can be performed at the frequency of the input signal or the frequency of the SSI mode. Similar to shallow foundations, charts or analytical equations can be developed to assist practitioners.

(iv) failure envelope identification: $a, b, \psi$ can be found using Eq. (9) and Table 1, while $c, d, e, f$ can be set to 1. They can be adjusted if experimental or numerical data are available.

(v) $c_1$: It mainly controls the horizontal-rotational behaviour and it is typically set to 1. It can be adjusted if experimental or numerical data are available.

(vi) $c_2, c_3$: These parameters primarily govern settlement and can typically be set to a value of 4.5. However, in scenarios involving a sequence of earthquake signals where the soil has achieved a sufficiently consolidated state—thereby limiting further settlement—a lower value, approximately 2.5, may be appropriate, as shown in Section 6.3, see also the work of Grange et al. [19] for a similar phenomenon during the TRISEE experimental campaign on shallow foundations. Adjustments can be made if experimental or numerical data are available.

(vii) $\kappa, \xi$: Set to 1 to ensure a horizontal tangent at $v = 0$ and to avoid counterintuitive vertical displacements. They can be adjusted if experimental or numerical data are available.

The simple and complex calibration methodologies for the macro-element parameters are summarized in Table 3. The coefficients $c_2$ and $c_3$ regulate the contribution of horizontal displacement $|u|$ and rotation $|\theta|$ to the vertical hardening variable $\gamma$, and therefore control the magnitude of cumulative settlement. Because settlement is highly sensitive to the cycle content and duration of the input record, these parameters cannot be uniquely identified across all six motions. A pragmatic procedure was adopted: stiffness and cyclic parameters were fixed by matching $M(\theta)$, after which $c_2, c_3$ were adjusted to reproduce the observed maximum settlement for each record. This record-conditioned calibration underlies Table 3 lists suggested values per case. A single global pair ($c_2, c_3$) can be used, but at the expense of a systematic settlement bias for some records. Consequently, the reported ranges should be interpreted as indicative rather than definitive, and further validation on an expanded dataset is needed before transferring these coefficients to other loading scenarios.

## 6. Validation

### 6.1. Dynamic centrifuge tests

To validate the proposed macro-element, results from the dynamic centrifuge experimental campaign conducted by Soriano et al. [38], Nohra et al. [39] and Nohra [54] are utilized. The primary objective of this campaign was to examine the dynamic response of a structure supported by a RI foundation. The structure in question was a cylinder with prototype-scale dimensions: a diameter of 6 m, a height of 12 m, and a total mass (including both the foundation and the structure) of 164.7 tons. It was positioned on a laminar container to minimize
<caption>**Table 3**</caption>

Calibration of the macro-element parameters (defaults in **bold**). Parenthetical values indicate record-conditioned alternatives identified from sequence-specific settlement response.

<table><thead><tr><th>Parameter</th><th>Simple (default)</th><th>Advanced (FE/impedance)</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>ASIRI [55] method: V<sub>0</sub> = A<sub>R</sub>q<sub>R</sub> + (1 - A<sub>R</sub>)q<sub>S</sub></td><td>3D FE, tangent method on monotonic V·w (bearing at first tangency)</td></tr><tr><td>k<sub>z</sub><sup>el</sup></td><td>Static FE: impose small u<sub>z</sub> on RI model, take slope F<sub>z</sub>/u<sub>z</sub> at origin</td><td>FE impedance on layered RI model: use ℜ{Z<sub>zz</sub>(ω<sub>SSI</sub>)} at the SSI frequency</td></tr><tr><td>k<sub>x</sub><sup>el</sup></td><td>Eq. (6) for circular footing,</td><td>.$$ \mathfrak{R}\{Z_{xx}(\omega_{SSI})\} $$</td></tr><tr><td>k<sub>θ</sub><sup>el</sup></td><td>Static FE: impose small θ, take slope M/θ at origin</td><td>.$$ \mathfrak{R}\{Z_{θθ}(\omega_{SSI})\} $$</td></tr><tr><td>C<sub>z</sub></td><td>Single dashpot at rep. SSI freq.: C<sub>z</sub> = ℑ{Z<sub>zz</sub>(ω<sub>SSI</sub>)}/ω<sub>SSI</sub></td><td>Multi-term ("monkey-tail") fit of Z<sub>zz</sub>(ω) over band of interest</td></tr><tr><td>C<sub>x</sub></td><td>C<sub>x</sub> = ℑ{Z<sub>xx</sub>(ω<sub>SSI</sub>)}/ω<sub>SSI</sub></td><td>As above; use same ω<sub>SSI</sub> as SSI mode</td></tr><tr><td>C<sub>θ</sub></td><td>C<sub>θ</sub> = ℑ{Z<sub>θθ</sub>(ω<sub>SSI</sub>)}/ω<sub>SSI</sub></td><td>As above</td></tr><tr><td>a, b, ψ</td><td>Eq. (9) with Table 1 (α, H<sub>LTP</sub>, φ<sub>LTP</sub>)</td><td>Fit to 3D FE monotonic sweeps (or tests) in (v, h, m)</td></tr><tr><td>c, d, e, f</td><td>1, 1, 1, 1</td><td>Adjust only if sweep data indicate different v-dependence</td></tr><tr><td>c<sub>1</sub></td><td>1</td><td>Adjust only if u<sup>pl</sup>-driven coupling to γ is mismatched</td></tr><tr><td>c<sub>2</sub>, c<sub>3</sub></td><td>4.5, 4.5 (2.5, 2.5 if prior consolidation limited settlement)</td><td>Adjust if cumulative settlement envelope change depending on record sequence</td></tr><tr><td>κ, ξ</td><td>1, 1</td><td>Adjust only to remove uplift artefacts in non-associated flow</td></tr></tbody></table>

Notes: (i) ω<sub>SSI</sub> denotes the predominant soil-foundation-structure interaction frequency (dominant or first coupled mode); (ii) values a, b, ψ in Eq. (9) are valid for circular RI-reinforced footings, a∈[2.5, 7]%, H<sub>LTP</sub>∈[0.5, 1.0] m, φ<sub>LTP</sub>∈[30°, 42°].

<caption>**Table 4**</caption>

Elastic and dynamic parameters of the centrifuge soil layers (after Alcala-Ochoa et al. [42]).

<table><thead><tr><th>Material/Layer</th><th>Model</th><th>ρ (kg/m³)</th><th>φ (°)</th><th>ψ (°)</th><th>V<sub>z</sub> (m/s)</th><th>E (MPa)</th></tr></thead><tbody><tr><td>LTP (Hostun sand mix)</td><td>Drucker-Prager</td><td>1600</td><td>36</td><td>2.5</td><td>0-72</td><td>7-26</td></tr><tr><td>Soft soil (80% kaolin + 20% sand)</td><td>Modified Cam-Clay</td><td>870</td><td>-</td><td>-</td><td>84-115</td><td>39-65</td></tr><tr><td>Stiff sand (HN31)</td><td>Drucker-Prager</td><td>1130</td><td>36</td><td>2.5</td><td>146-160</td><td>93-117</td></tr></tbody></table>

boundary effects. The coverage area α was 7% and 7 RI were employed. A three-layer soil profile was selected for the experiment, consisting of, from top to bottom: dry Hostun sand mix for the load transfer platform (LTP), Speswhite kaolin mix, and saturated Hostun sand, Fig. 7.

The centrifuge soil model reproduces the ASIRIplus-SDS configuration described in the work of Nohra et al. [54,56] for dynamic tests. It consists of three layers: (i) 1 m Load Transfer Platform (LTP) of dense Hostun sand (mix of HN38, HN34, HN31, HN04/08 and HN06/1), (ii) 9.5 m thick soft soil composed of 80% Speswhite kaolin and 20% Fontainebleau sand, and (iii) 4 m stiff bearing layer of dense Hostun sand (HN31). The soft layer was pre-consolidated at 120 kPa to achieve an undrained shear strength of about 20 kPa at 1 g. The water table is located at the base of the LTP, 1 m below the footing level. The principal geotechnical and dynamic parameters adopted in the numerical simulations are summarized in Table 4. These values are consistent with those used in the work of Alcala-Ochoa et al. [42], which were calibrated against laboratory and centrifuge data from the ASIRIplus-SDS campaign.

To assess the soil response, accelerometers (and pore pressure transducers) were positioned within the laminar container. The structure was also equipped with accelerometers in both horizontal and vertical directions to monitor the response of the rigid cylinder structure, and laser sensors to track displacements. The positions of the sensors are shown in Fig. 8. For simplicity, only the sensors related to the validation of the macro-element are labelled in the figure.

Accelerometer A1 records the input signal at the base of the container. The overturning moment at the base of the foundation was calculated using accelerometer A5 (by multiplying the recorded acceleration with the mass and height of the structure). The rotation at the base was calculated from A5 and A7, after double integration to obtain displacements D5 and D7, followed by calculating their difference and dividing by the height of the structure. Settlement was calculated as the average of the vertical displacements D1 and D2. The acceleration recorded by A19, positioned near the soil surface and at a considerable distance from the structure, is considered the ground surface acceleration and the input signal to the macro-element.

The container with the soil and structure was placed in the centrifuge, and the gravitational force was gradually increased from 0 to 50 g. After nearly two hours of consolidation in flight, the model

<caption>**Table 5**</caption>

Signal type and PGA of the sequence's input motions [39,54].

<table><thead><tr><th>Signal</th><th>PGA (g)</th></tr></thead><tbody><tr><td>Landers 0.05 g</td><td>0.72</td></tr><tr><td>Northridge 0.05 g</td><td>0.074</td></tr><tr><td>Sinus 50 Hz 0.05 g</td><td>0.119</td></tr><tr><td>Sinus 50 Hz 0.15 g</td><td>0.261</td></tr><tr><td>Sinus 50 Hz 0.25 g</td><td>0.364</td></tr><tr><td>Landers 0.05 g</td><td>0.079</td></tr><tr><td>Sinus 90 Hz 0.05 g</td><td>0.086</td></tr><tr><td>Sinus 90 Hz 0.15 g</td><td>0.275</td></tr><tr><td>Sinus 90 Hz 0.25 g</td><td>0.397</td></tr><tr><td>Landers 0.05 g</td><td>0.075</td></tr><tr><td>Sinus 120 Hz 0.05 g</td><td>0.141</td></tr><tr><td>Sinus 120 Hz 0.15 g</td><td>0.292</td></tr><tr><td>Sinus 120 Hz 0.25 g</td><td>0.412</td></tr><tr><td>Landers 0.05 g</td><td>0.08</td></tr><tr><td>Landers 0.15 g</td><td>0.196</td></tr><tr><td>Landers 0.30 g</td><td>0.328</td></tr><tr><td>Northridge 0.05 g</td><td>0.084</td></tr><tr><td>Northridge 0.15 g</td><td>0.225</td></tr><tr><td>Northridge 0.30 g</td><td>0.464</td></tr><tr><td>Landers 0.05 g</td><td>0.075</td></tr></tbody></table>

stabilized. Then, twenty ground motions were applied at the base of the model, as detailed in Table 5.

Typical experimental results are illustrated in Fig. 9, as referenced in the report of Nohra et al. [39]. These findings demonstrate an increase in acceleration at the top of the structure, the progression of settlement, and the dissipation and nonlinearities observed in the moment-rotation curve. The results presented pertain specifically to the Landers earthquake with a peak ground acceleration of 0.15 g.

Although the centrifuge campaign included twenty successive ground motions, each numerical simulation in the present study was performed independently, without considering cumulative effects. This assumption is supported by the experimental observations of Nohra [54], who reported only a moderate and largely reversible degradation of the small-strain stiffness after each excitation. Measurements from bender elements showed that the maximum shear modulus (G<sub>max</sub>) experienced their largest decrease during the first strong motions (1 Hz sine series), followed by partial recovery and stabilization for the subsequent earthquakes.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Schematic representation of the centrifuge model, Nohra *et al.* [39,54].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Centrifuge model of the dynamic test of a RI foundation, prototype scale: plane and front view, adopted from Nohra *et al.* [39,54].</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Experimental results, Landers 0.15 g: acceleration (a) at the ground surface, (b) at the top of the structure, (c) settlement vs. moment and (d) moment vs. rotation [39,54].</figcaption></figure>

The total reduction of $G_{max}$ throughout the entire sequence did not exceed about 30%, indicating that the soil stiffness at small strains was almost recovered between successive motions. Furthermore, the evolution of excess pore pressures in the clay layer revealed limited accumulation and a progressive stabilization between events, suggesting that cyclic degradation and pore-pressure buildup were not significant over the complete series. Consequently, the independent treatment of each input motion in the numerical analyses is deemed consistent with the experimentally observed behaviour of the reinforced soil system.

## 6.2. Numerical model

The numerical model consists of a simplified system in which the RI foundation and the surrounding soil are represented by a single ME positioned at the base of the structure, as shown in Fig. 10. The superstructure is modelled as a beam with a lumped mass applied at its top. The six earthquake signals highlighted in bold in Table 5 were applied independently, with no cumulative effects considered, as justified in Section 6.1. Numerical simulations were performed in the finite element code Cast3M [52].

The exact position of the centre of mass of the rigid cylinder was estimated based on the dynamic centrifuge tests [39,54]. In the numerical model, the lumped mass was placed at $H = 9$ m, which ensures that the first numerical eigenfrequency matches the experimentally identified fundamental frequency of 8.0 Hz [54]. This procedure follows standard calibration practices in substructured soil-structure interaction analyses, where the elevation of the concentrated mass is adjusted to reproduce the measured modal properties without modifying the foundation stiffness parameters. Within the experimental uncertainty, this height reproduces the observed modal frequency and shape while preserving the overall stiffness characteristics and global dynamic behaviour of the system.

The finite element code Cast3M [52] was used for the numerical calculations. The model was solved using the implicit Newmark scheme with a constant time step $\Delta t = 0.0052$ s. A Rayleigh damping ratio of 2% was applied to the structure, while hysteretic energy dissipation arose naturally from the macro-element plasticity. Nonlinear equilibrium was achieved through Newton-Raphson iterations with the consistent tangent stiffness, using Cast3M's default convergence tolerances (relative residuals below $10^{-6}$). The system was implemented as a zero-length

<table><caption>**Table 6**
Parameters of the macro-element.</caption><thead><tr><th>Parameter</th><th>Value</th><th>SI Unit</th></tr></thead><tbody><tr><td>V<sub>0</sub></td><td>4593 × 10<sup>3</sup></td><td>N</td></tr><tr><td>k<sub>z</sub><sup>el</sup></td><td>3.98 × 10<sup>8</sup></td><td>N m<sup>-1</sup></td></tr><tr><td>k<sub>x</sub><sup>el</sup></td><td>1.45 × 10<sup>8</sup></td><td>N m<sup>-1</sup></td></tr><tr><td>k<sub>θ</sub><sup>el</sup></td><td>1.4 × 10<sup>9</sup></td><td>N rad<sup>-1</sup></td></tr><tr><td>a</td><td>0.7064</td><td>–</td></tr><tr><td>b</td><td>0.442</td><td>–</td></tr><tr><td>c, d, e, f</td><td>1 1 1 1</td><td>–</td></tr><tr><td>ψ</td><td>-0.3325</td><td>–</td></tr><tr><td>c<sub>1</sub></td><td>1</td><td>–</td></tr><tr><td>c<sub>2</sub>, c<sub>3</sub> (Landers)</td><td>4.5 4.5</td><td>–</td></tr><tr><td>c<sub>2</sub>, c<sub>3</sub> (Northridge)</td><td>2.5 2.5</td><td>–</td></tr><tr><td>κ, ξ</td><td>1 1</td><td>–</td></tr></tbody></table>

joint element (JOI1) representing the macro-element and a two-node beam element (POUT) representing the structure. The base node was fully fixed, and the upper node carried the lumped mass and horizontal inertial loading derived from the experimental accelerations. The macro-element's parameters are summarized in Table 6.

The dynamic analysis was performed using the PASAPAS procedure, with forces, moments, displacements, and rotations recorded at both nodes for post-processing and comparison with experimental data.

* $V_0$ is calculated using a 3D FEM model and the tangent method [41,42].

* Impedance calculations were performed using a 3D multi-layer finite FEM model, following the procedure detailed in the work of Coronado and Gidwani [57]. The results for the dynamic stiffnesses are presented in Fig. 11. The stiffness terms $k_z^{el}$, $k_x^{el}$, and $k_θ^{el}$ were extracted at the system SSI frequency (approximately 1 Hz), which lay within a low-frequency plateau of $\Re\{Z(\omega)\}$ (0–2 Hz).

* Radiation damping is not taken into account in this analysis, as the dynamic centrifuge tests ensured that all stress wave energy was contained within the laminar container.

* $a$, $b$, and $\psi$ are calculated using Eq. (9) and Table 1, considering $\alpha = 7\%$, $H_{LTP} = 1$ m and $\varphi_{LTP} = 36^\circ$.

* For all the Landers signals $c_1 = 1$ and $c_2 = c_3 = 4.5$.

* For all the Northridge signals $c_1 = 1$ and $c_2 = c_3 = 2.5$, due to consolidation (see discussion in Section 5).

The parameter $V_0$ was determined from three-dimensional finite element calculations using the tangent method [41,42]. The elastic stiffnesses $k_z^{el}$, $k_x^{el}$, and $k_θ^{el}$ were extracted from dynamic impedance analyses of a multilayer model. The impedance results are shown in Fig. 11, where the stiffness values at zero frequency were adopted. Radiation damping was neglected, as the dynamic centrifuge experiments ensured that stress waves were confined within the laminar container. The parameters $a$, $b$, and $\psi$ were evaluated using Eq. (9) and Table 1, for $\alpha = 7\%$, $H_{LTP} = 1$ m, and $\varphi_{LTP} = 36^\circ$. The coefficients $c_1$, $c_2$, and $c_3$ were selected following the calibration procedure in Section 5: $c_1 = 1$, and $(c_2, c_3) = (4.5, 4.5)$ for the Landers records and $(2.5, 2.5)$ for the Northridge records, the latter reflecting partial consolidation between successive excitations.

## 6.3. Numerical vs. experimental results

The ME predictions are compared with centrifuge data [39] for six input motions: Landers (0.05 g, 0.15 g, 0.30 g) and Northridge (0.05 g, 0.15 g, 0.30 g). Each comparison figure contains three panels: (a) relative top displacement time history $u_{top}(t)$, (b) settlement-rotation $\omega(\theta)$, and (c) moment-rotation $M(\theta)$. Legends are unified to Experimental
<figure><img src="image_13.png" /><figcaption>Fig. 10. Numerical model for dynamic analyses: (a) RI foundation represented by a ME combining nonlinear behaviour, dynamic stiffness, and radiation damping; (b) equivalent lumped-mass system subjected to seismic input.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Frequency-dependent dynamic stiffnesses for horizontal, vertical, and rotational degrees of freedom, obtained from FEM impedance calculations.</figcaption></figure>

(grey) and Macro-element (blue). To complement the qualitative over-lays, two quantitative metrics are reported in each subplot: the squared correlation $R^2$ (shape agreement, $0 \le R^2 \le 1$) and the Normalized Mean Absolute Error NMAE (%), defined as:

$$NMAE = 100 \frac{\frac{1}{n} \sum_i |y_i^{\text{exp}} - y_i^{\text{num}}|}{\max(y^{\text{exp}}) - \min(y^{\text{exp}})} \quad (30)$$

with normalization by the range of the experimental ordinate so that values are comparable across panels and cases (smaller percentages indicate better agreement).

The relative top displacement $u_{\text{top}}(t)$ shows moderate agreement between experimental and numerical results, with $R^2$ values ranging from about 0.37 to 0.64 and NMAE values between 4.5% and 9%. These levels of discrepancy are consistent with a residual mismatch in the effective soil-structure interaction frequency and damping characteristics, which slightly affect phase and amplitude in the simulated time histories.

The moment-rotation response $M(\theta)$ displays the best overall agreement. In four of the six cases, $R^2$ ranges between approximately 0.76 and 0.91 with modest amplitude error NMAE between 5% and 11%. The cyclic shape and hysteresis loops are reproduced satisfactorily, with only small constant biases of the order of $\pm(0.01-0.06)$ MN m depending on the record. This behaviour indicates that the macro-element reproduces the rotational stiffness and cyclic degradation of

the foundation-soil system with reasonable fidelity at the loading levels tested.

For the settlement-rotation behaviour $\omega(\theta)$, the metric $R^2$ was deliberately omitted because it is not informative in this context. The $\omega(\theta)$ plots form multi-valued, looped trajectories parameterized by time, whereas $R^2$ assumes a single-valued functional relation and linear correlation of paired samples. Small phase shifts during load reversals or accumulation of permanent offsets can therefore collapse $R^2$ to near zero even when the overall envelopes are comparable. In panel (b), only the NMAE and qualitative agreement of amplitude and envelope are considered.

A consistent underestimation of settlement magnitude is observed. The corresponding NMAE values are comparatively large, about 25%–36%, with the exception of Northridge 0.30 g where NMAE $\approx 4.9\%$. This bias suggests that the vertical plastic response remains too stiff and that the coupling between vertical and horizontal-rotational plastic mechanisms in the current formulation is insufficient.

Taken together, these results indicate that the macro-element captures the essential features of the nonlinear dynamic response: the rotational stiffness and cyclic shape are reproduced accurately, while the horizontal and vertical components still show moderate deviations. The moderate accuracy of $u_{\text{top}}(t)$ is consistent with minor differences. The underestimation of settlement is attributed to the current evolution law of the variable $\gamma$, which governs vertical hardening. Targeted adjustment of the coefficients governing the expansion of the loading surface in the vertical direction (parameters controlling the contribution of $|u|$, $|\theta|$ to the evolution of $\gamma$) is expected to reduce this bias, while leaving the moment-rotation response largely unaffected. In practice, the coefficients $c_2$ and $c_3$ were overhauled to reproduce the observed maximum settlement in each record. Because vertical accumulation depends strongly on record-specific cycle content and duration, equal-PGA motions such as Landers and Northridge yield different optimal ($c_2, c_3$). For this reason, Table 3 reports recommended, record-consistent values rather than unique optima from a global identification. These values should be regarded as case-specific calibrations within the tested motion set (see Figs. 12–17); further studies with a broader suite of records are required before general design recommendations can be drawn.

### 6.4. Frequency-domain comparison and floor spectra

The frequency-domain comparison relies on: (i) the single-sided FFT amplitude of the absolute top displacement, $|\hat{u}_{\text{top}}(f)|$, and (ii) the 2%-damped floor spectrum $S_d(f)$. A light damping level ($\xi = 2\%$)
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Validation against centrifuge data [39] for Landers 0.05 g. Panels: (a) relative top displacement $u_{top}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1$, $c_2 = c_3 = 4.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Validation against centrifuge data [39] for Landers 0.15 g. Panels: (a) relative top displacement $u_{top}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1$, $c_2 = c_3 = 4.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

is adopted so that $S_d(f)$ retains spectral detail-minimizing smoothing and preserving narrow-band SSI features. Six two-panel figures are provided (Experimental vs. Macro-element), one per input motion (Landers 0.05 g, 0.15 g, 0.30 g; Northridge 0.05 g, 0.15 g, 0.30 g). Panel (a) shows $|\hat{u}_{top}(f)|$ normalized by the case-specific SSI peak amplitude of the experimental response, so that spectral shape and peak frequency are emphasized rather than absolute height. Panel (b) reports $S_d(f)$ for the top absolute motion in meters, with the surface input spectrum for

reference. Unless noted otherwise, amplitude statements below refer to panel (b), $S_d$.

The laminar container used in the centrifuge campaign confines wave energy, radiation damping is therefore excluded from the time-domain model. Dissipation then arises from the ME plasticity and 2% Rayleigh damping in the superstructure, with stiffness terms taken from the impedance-based calibration introduced earlier around the SSI band. Under these assumptions the response is governed by the coupled
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Validation against centrifuge data [39] for Landers 0.3 g. Panels: (a) relative top displacement $u_{\text{top}}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1, c_2 = c_3 = 4.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Validation against centrifuge data [39] for Northridge 0.05 g. Panels: (a) relative top displacement $u_{\text{top}}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1, c_2 = c_3 = 2.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

foundation-superstructure mode and by the spectral content of each record.

For Landers, experimental floor spectra peaks lie near 0.67–0.75 Hz; for Northridge, peaks group around 0.83–1.03 Hz. This difference comes directly from the input content that each record supplies to the SSI band and from the effective coupled frequency: Landers carries stronger long-period energy and preferentially excites the lower part of the band,

while Northridge contains more energy near ~ 1 Hz and activates a slightly stiffer branch of the coupled response. This reading is consistent with the time-domain comparisons (Section 6.3) where Landers motions drive larger quasi-rigid components and Northridge excites a tighter oscillatory cycle.

Across all six records, the top floor spectrum exceeds the surface input spectrum when evaluated at the respective peaks: ≈ 1.84
<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Validation against centrifuge data [39] for **Northridge 0.15 g.** Panels: (a) relative top displacement $u_{\text{top}}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1$, $c_2 = c_3 = 2.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Validation against centrifuge data [39] for **Northridge 0.3 g.** Panels: (a) relative top displacement $u_{\text{top}}(t)$, (b) settlement-rotation $\omega(\theta)$, (c) moment-rotation $M(\theta)$ ($c_1 = 1$, $c_2 = c_3 = 2.5$). Experimental (grey) and macro-element (blue). The legend reports $R^2$ (shape agreement) and NMAE (%) referenced to the experimental range.</figcaption></figure>

(Landers 0.05 g), 2.94 (Landers 0.15 g), 2.54 (Landers 0.30 g), 3.44 (Northridge 0.05 g), 2.98 (Northridge 0.15 g), and 2.66 (Northridge 0.30 g). These values quantify SSI-induced amplification: the top motion behaves as a selective filter that concentrates energy around the coupled mode once rigid-body and rocking components are transmitted to the structure.

Peak locations and amplitudes at the top indicate the following patterns:

(i) Landers 0.05 g. Experimental peak at 0.75 Hz (0.083 m); numerical peak at 1.36 Hz (0.048 m). The upward shift (+81%) indicates an effective stiffness that is too high in this low-amplitude regime; the amplitude deficit points to excess damping or insufficient hysteretic mobilization at small strains.

(ii) Landers 0.15 g and 0.30 g. Peak frequencies align (~0.67–0.69 Hz), meaning the macro-element reproduces the coupled period in the Landers band. The lower numerical amplitudes (ratios ≈
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Landers 0.05 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{\text{top}}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. Frequency axis 0–5 Hz; common time base. The ME peak frequency exceeds the experimental one, and the top $S_d$ amplitude is lower (panel (b)).</figcaption></figure>

0.38 and 0.69) indicate stronger effective dissipation or reduced inelastic drift compared with the experiment; this is consistent with the time-domain results.

(iii) Northridge 0.05 g–0.30 g. Numerical peaks occur at slightly higher frequencies (+3.5%, +6.4%, +22% for 0.05 g, 0.15 g, 0.30 g) and exceed experimental amplitudes (ratios ≈ 1.24, 1.14, 1.33). The combination of a modestly stiffer effective response and lower hysteretic loss in the model explains the upward drift in $f_{\text{peak}}$ and the amplitude surplus.

The top-over-surface amplification and the narrow peaks indicate resonance in the SSI-relevant band. With radiation damping absent in the laminar container, wave energy remains within the soil-foundation-structure system. In the numerical model, peak height and bandwidth are then governed by the structural damping and ME hysteresis. For Landers (0.15 g–0.30 g), the ME matches the peak frequency while the top $S_d$ amplitudes remain lower than in the experiment, consistent with a more dissipative effective response in the model. For Northridge (0.05 g–0.30 g), the ME peaks occur at slightly higher frequency and the $S_d$ amplitudes exceed the experiment, consistent with a stiffer and less dissipative effective response around ~1 Hz.

Low-amplitude shifts (e.g., Landers 0.05 g) are particularly sensitive to the small-strain elastic stiffness embedded in the macro-element and to the structural mass-height idealization used to match the modal properties; they are also affected by the relative weighting of hysteretic and structural damping around the identified SSI band. These sensitivities already appear in the time-domain comparisons and remain consistent with the calibration pathway adopted for the laminar container configuration (see Figs. 18–23).

## 7. Conclusions

A plasticity-based macro-element (ME) was developed for circular shallow foundations on rigid-inclusion (RI)-reinforced, layered ground and was validated under seismic loading conditions. The formulation combined: (i) an inclined interaction ellipse in the $(h, m)$ plane, sized and tilted through $(a, b, \psi)$ linked to RI descriptors $(\alpha, H_{\text{LTP}}, \varphi_{\text{LTP}})$ within a stated validity domain; (ii) elastic terms extracted from impedance functions at the soil-structure interaction frequency $\omega_{\text{SSI}}$; and (iii) a vertical hardening variable $\gamma$ whose evolution included explicit contributions from plastic horizontal displacement and rotation. A diagonal elastic matrix was adopted as standard for shallow circular footings with small eccentricities. This first-order approximation was consistent with the axial symmetry of the configuration, and it enabled a practical calibration without degrading the match to the centrifuge observations for the cases analysed.

The ME reproduced key response features observed in dynamic centrifuge tests on RI-reinforced ground. Across six records (Landers and Northridge at three PGAs each): (i) the moment-rotation loops $M(\theta)$ were captured with good fidelity; (ii) relative top displacement time histories showed moderate agreement; and (iii) cumulative settlement was generally underestimated, indicating that the current $\gamma$-evolution limited vertical accumulation. Within the stated domain, the ME was an efficient tool for nonlinear SSI analyses of RI-supported structures, enabling performance-oriented checks with modest calibration effort.

A pragmatic calibration was established: $V_0$ from ASIRI [55] (with optional FEM confirmation); $k_x^{\text{el}}$ from the Gazetas closed form for circular footings on homogeneous profiles (or from impedance at $\omega_{\text{SSI}}$ when layering or RI effects are significant); $(k_z^{\text{el}}, k_\theta^{\text{el}})$ from impedance functions at $\omega_{\text{SSI}}$; radiation terms $C_{ii}$ likewise from impedance at $\omega_{\text{SSI}}$ for field applications, whereas in the present laminar-container centrifuge configuration far-field radiation was negligible and energy dissipation was adequately represented by macro-element hysteresis and structural damping, so explicit radiation dashpots were not included; $(a, b, \psi)$ from the closed-form fit when within its domain; and record-conditioned reductions of $(c_2, c_3)$ when prior consolidation limited further settlement.

The inclined $(h, m)$ envelope and its calibration were limited to circular footings and to the ranges used to build the fit: $\alpha \in [2.5, 7.0]\%$, $H_{\text{LTP}} \in [0.5, 1.0]$ m, and $\varphi_{\text{LTP}} \in [30^\circ, 42^\circ]$. Use outside these ranges requires new calibration (numerical sweeps or tests). The elastic matrix was taken diagonal; inclusion of coupling may be required for embedded or skirted bases, markedly non-circular plan shapes, large eccentricities, strong layering/anisotropy or sloping ground, or when the response is governed by high-frequency rocking near an SSI mode. Settlement predictions were most sensitive to the $\gamma$-evolution; the present form provided control via $(c_2, c_3)$ but remained a first-order approximation.

Compared with classical shallow-footing macro-elements that used a *rugby-ball* ultimate surface and case-by-case calibration, and with hypoplastic variants for sands, the present RI-ME introduced three practical advances: (i) a closed-form inclined $(h, m)$ envelope tied to RI descriptors within stated limits, enabling direct use without a sweep campaign; (ii) linear terms taken from impedances at $\omega_{\text{SSI}}$ rather than zero-frequency stiffness or frequency-independent dashpots; and (iii) calibration/validation against dynamic centrifuge tests on RI-reinforced, layered ground.

Four developments are recommended: (i) extend the envelope calibration to broader $(\alpha, H_{\text{LTP}}, \varphi_{\text{LTP}})$ ranges and to non-circular foundations; (ii) refine the $\gamma$-evolution to improve settlement accumulation-including cycle-content sensitivity-and, where needed, introduce directional coupling via off-diagonal elastic terms inferred from impedance; (iii) incorporate strain-dependent material damping, for example
<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Landers 0.15 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{top}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. The ME reproduces the experimental peak frequency; the top $S_d$ amplitude is lower, panel (b).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Landers 0.30 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{top}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. The ME matches the experimental peak frequency; the top $S_d$ amplitude is lower, panel (b).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Northridge 0.05 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{top}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. The ME peak occurs at a slightly higher frequency and the top $S_d$ amplitude is higher than the experiment, panel (b).</figcaption></figure>

through a simple hyperbolic stress–strain relation or an equivalent damping model; and (iv) compile charts or rational fits for $K_{ii}(\omega)$ and

$C_{ii}(\omega)$ for representative RI-layered profiles to reduce the need for on project-specific impedance analyses.
<figure><img src="image_25.png" /><figcaption>**Fig. 22.** Northridge 0.15 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{top}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. The ME peak frequency is slightly higher and the top $S_d$ amplitude is higher, panel (b).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Northridge 0.30 g-frequency-domain comparison. (a) Normalized $|\hat{u}_{top}(f)|$. (b) $S_d(f)$ and surface input spectrum. Frequency axis 0–5 Hz. The ME peak shifts further to higher frequency and the top $S_d$ amplitude is higher than the experiment (panel (b)).</figcaption></figure>

### CRediT authorship contribution statement

**Ramon Alcala-Ochoa:** Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation. **Zheng Li:** Writing – review & editing, Supervision, Investigation, Conceptualization. **Panagiotis Kotronis:** Writing – review & editing, Supervision, Investigation, Funding acquisition, Conceptualization. **Giulio Sciarra:** Writing – review & editing, Supervision, Investigation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

The financial support from the French National Research Agency for the national project: ASIRIplus-SDS (Amélioration des sols par inclusoins rigides : Sollicitations Dynamiques et Sismiques) with NO. Grant ANR-19-CE22-0015 is greatly acknowledged. The authors would like to acknowledge the help from the SuperComputing Center of École Centrale de Nantes.

### Data availability

Data may be available upon request.

### References

[1] Mayoral J, Romo M, Cirion A, Paulin J. Effect of layered clay deposits on the seismic behavior of a rigid inclusion. In: Proceedings of the symposium on rigid inclusions in difficult subsoil condition. Vol. 1, 2006, p. 1–12.

[2] Hatem A. Comportement en zone sismique des inclusions rigides : Analyse de l'interaction sol-inclusion-matelas de répartition - structure. (Ph.D. thesis), University Lille 1; 2009, http://www.theses.fr/2009LIL10097.

[3] Diaz-Guzman M, Manica-Malcom M, Ovando S, et al. Respuesta sísmica de sistemas de cimentación a base de inclusiones rígidas. In: XXX Reunión Nacional de Ingeniería Geotécnica Respuesta sísmica de sistemas de cimentación a base de inclusiones rígidas. 2021.

[4] López-Jiménez GA. Static and dynamic behaviour of pile supported structures in soft soil. (Ph.D. thesis), Université Grenoble Alpes; 2019, https://theses.hal.sciences-tel-02271948.

[5] Saadatinezhad M, Lakirouhani A, Asli SJabini. Seismic response of non-connected piled raft foundations. Int J Geotech Eng 2021;15:66–80. http://dx.doi.org/10.1080/19386362.2019.1565392.

[6] Mánica Malcom MÁ, Ovando-Shelley E, Botero Jaramillo E. Numerical study of the seismic behavior of rigid inclusions in soft Mexico city clay. J Earthq Eng 2016;20:447–75. http://dx.doi.org/10.1080/13632469.2015.1085462.

[7] López Jiménez GA, Dias D, Jenck O. Effect of layered liquefiable deposits on the seismic response of soil-foundations-structure systems. Soil Dyn Earthq Eng 2019;124:1–15. http://dx.doi.org/10.1016/j.soildyn.2019.05.026, URL: https://www.sciencedirect.com/science/article/pii/S0267726118313654.

[8] J.L. Rangel-Núñez, Gómez-Bernal A, Aguirre-González J, Sordo-Zabay E, Ibarra-Razo E. Dynamic response of soft soil deposits improved with rigid inclusions. In: 14th world conference on earthquake engineering. 14WCEE, 2008, p. 1–8.

[9] Han X, Li Y, Ji J, Ying J, Li W, Dai B. Numerical simulation on the seismic absorption effect of the cushion in rigid-pile composite foundation. Earthq Eng Vib Vib 2016;15:369–78. http://dx.doi.org/10.1007/s11803-016-0329-x.
[10] Hu Z, Li J, Lin G. Seismic nonlinear interaction and uplifting effects between the nuclear island building and soil. Eng Struct 2025;322:119168. http://dx.doi.org/10.1016/j.engstruct.2024.119168.

[11] Perez-Herreros J. Dynamic soil-structure interaction of pile foundations: experimental and numerical study. (Ph.D. thesis), École Centrale de Nantes.; 2020. https://tel.archives-ouvertes.fr/tel-02972312.

[12] Nova R, Montrasio L. Settlements of shallow foundations on sand. Géothéchnique 1991;41:243-56.

[13] Cremer C. Modélisation du comportement non linéaire des fondations superficielles sous séisme: macro-élément d'interaction sol-structure. (Ph.D. thesis), Cachan, École normale supérieure; 2001.

[14] Cremer C, Pecker A, Davenne L. Cyclic macro-element for soil-structure interaction: material and geometrical non-linearities. Int J Numer Anal Methods Geomech 2001;25:1257-84. http://dx.doi.org/10.1002/NAG.175.

[15] Cremer C, Pecker A, Davenne L. Modelling of nonlinear dynamic behaviour of a shallow strip foundation with macro-element. J Earthq Eng 2002;6:175-211. http://dx.doi.org/10.1080/13632460209350414.

[16] Chatzigogos C, Pecker A, Salençon J. A macro-element for dynamic soil-structure interaction analyses of shallow foundations. In: 4th international conference of earthquake and geotechnical engineering. Greece: Thessaloniki; 2007. URL: https://hal.science/hal-00263677. paper n° 1387.

[17] Grange S. Modélisation simplifiée 3D de l'interaction sol-structure: application au génie parasismique (Theses), Institut National Polytechnique de Grenoble - INPG; 2008. https://theses.hal.science/tel-00306842.

[18] Grange S, Kotronis P, Mazars J. A macro-element for a circular foundation to simulate 3d soil-structure interaction. Int J Numer Anal Methods Geomech 2008;32:1205-27. http://dx.doi.org/10.1002/nag.664.

[19] Grange S, Kotronis P, Mazars J. A macro-element to simulate 3d soil-structure interaction considering plasticity and uplift. Int J Solids Struct 2009;46:3651-63. http://dx.doi.org/10.1016/j.ijsolstr.2009.06.015.

[20] Grange S, Kotronis P, Mazars J. A macro-element to simulate dynamic soil-structure interaction. Eng Struct 2009;31:3034-46.

[21] Chatzigogos CT, Pecker A, Tamagnini C. A hypoplastic macroelement model for shallow foundations. Soil Dyn Earthq Eng 2009;29:765-81. http://dx.doi.org/10.1016/j.soildyn.2008.08.009.

[22] Salciarini D, Tamagnini C. A hypoplastic macroelement model for shallow foundations under monotonic and cyclic loads. Acta Geotech 2009;4:163-76. http://dx.doi.org/10.1007/s11440-009-0087-2.

[23] Salciarini D, Tamagnini C, Grange S, Kotronis P. La modellazione dei fenomeni di interazione 719terreno-struttura mediante macroelementi: elastoplasticità vs. ipoplasticità. Rivista Italiana di Geotech 2010;4:9-28.

[24] Wang Y. A macro-element model for spudcan foundations in clay overlying sand. (Ph.D. thesis), University of Western Australia.; 2021. https://research-repository.uwa.edu.au/en/publications/a-macro-element-model-for-spudcan-foundations-in-clay-overlying-s.

[25] da Silva AP, Diambra A, Karamitros D. Macro-element modelling of suction-embedded plate anchors for floating offshore structures. In: E3S web of conferences. Vol. 92, 2019, http://dx.doi.org/10.1051/E3SCONF/20199216009.

[26] Jin Z, Yin ZY, Kotronis P, Li Z. Advanced numerical modelling of caisson foundations in sand to investigate the failure envelope in the H-M-V space. Ocean Eng 2019;190). http://dx.doi.org/10.1016/j.oceaneng.2019.106394.

[27] Jin Z, Li Z, Yin ZY, Kotronis P. Numerical modeling of soil-pipe interaction of single pipeline at shallow embedment in clay by hypoplastic macroelement. Ocean Eng 2021;110017. http://dx.doi.org/10.1016/j.oceaneng.2021.110017.

[28] Correia A. A pile-head macro-element approach to seismic design of monoshaft-supported bridges. (Ph.D. thesis), ROSE Programme, UME School, IUSS Pavia; 2011.

[29] Correia A, Pecker A, Kramer S, Pinho R. Nonlinear pile-head macro-element model: Ssi effects on the seismic response of a monoshaft-supported bridge. In: Proceedings, 15th world conference on earthquake engineering. Lisbon, Portugal; 2012.

[30] Correia AA, Pecker A. Nonlinear pile-head macro-element for the seismic analysis of structures on flexible piles. Bull Earthq Eng 2021;19:1815-49. http://dx.doi.org/10.1007/s10518-020-01034-4.

[31] Gorini DN, Callisto L. A multiaxial inertial macroelement for deep foundations. Comput Geotech 2023;155. http://dx.doi.org/10.1016/j.compgeo.2022.105222.

[32] Li Z, Kotronis P, Escoffier S, Tamagnini C. A hypoplastic macroelement for single vertical piles in sand subject to three-dimensional loading conditions. Acta Geotech 2015;11. http://dx.doi.org/10.1007/s11440-015-0415-7.

[33] Li Z, Kotronis P, Escoffier S, Tamagnini C. A hypoplastic macroelement formulation for single batter piles in sand. Int J Numer Anal Methods Geomech 2018;42:1346-65. http://dx.doi.org/10.1002/nag.2794.

[34] Doulgeroglou AA, Kotronis P, Sciarra G, Bouillon C. 3D interaction diagrams for symmetrically reinforced concrete square sections with various reinforcement ratios. Eng Struct 2022;262:114272. http://dx.doi.org/10.1016/j.engstruct.2022.114272.

[35] Hill R. The mathematical theory of plasticity. Clarendon Press; 1950.

[36] Borja RI. Plasticity. Springer Berlin Heidelberg; 2013, doi 10 978-3.

[37] Grange S, Kotronis P, Mazars J. A macro-element for a shallow foundation to simulate soil-structure interaction considering uplift. Comptes Rendus Mécanique 2008;336:856-62.

[38] Soriano C, Escoffier S, Li Z. Data Report (2) ASIRI+ Project Amélioration Des Sols Par Inclusions Rigides: Sollicitations Dynamiques Et Sismiques. Technical Report., Université Gustave Eiffel; 2022.

[39] Nohra C, Escoffier S, Li Z, Soriano C. Data report (9) soil-structure interaction: building over rigid inclusions ASIRI+ project amélioration des sols par inclusions rigides: sollicitations dynamiques et sismiques. Technical Report., Université Gustave Eiffel; 2023.

[40] Gazetas G. Analysis of machine foundation vibrations: state of the art. Int J Soil Dyn Earthq Eng 1983;2:2-42.

[41] Alcala-Ochoa R. Numerical modelling of rigid inclusions foundations subjected to seismic loading (Ph.D. thesis), École centrale de Nantes; 2023, https://theses.hal.science/tel-04516907.

[42] Alcala-Ochoa R, Li Z, Kotronis P, Sciarra G. 3D failure envelope of rigid inclusion reinforced foundations. Acta Geotech 2024;19:6973-92. http://dx.doi.org/10.1007/s11440-024-02309-9.

[43] Gottardi G, Butterfield R. On the bearing capacity of surface footings on sand under general planar loads. Soils Found 1993;33:68-79. http://dx.doi.org/10.3208/sandf1972.33.3.68.

[44] Cassidy M. Non-linear analysis of jack-up structures subjected to random waves (Ph.D. thesis), University of Oxford.; 1999, https://ora.ox.ac.uk/objects/uuid:02b96107-638b-4200-9fd7-fa49635594e3.

[45] Gottardi G, Houlsby GT, Butterfield R. Plastic response of circular footings on sand under general planar loading. 1999, http://dx.doi.org/10.1680/GEOT.1999.49.4.453.

[46] Prevost J. Anisotropic undrained stress-strain behavior of clays. ASCE J Geotech Eng Div 1978;104:1075-90.

[47] Wolf JP. Soil-structure-interaction analysis in time domain. Nucl Eng Des 1989;111:381-93.

[48] Gazetas G. Formulas and charts for impedances of surface and embedded foundations. J Geotech Eng 1991;117:1363-81. http://dx.doi.org/10.1061/(ASCE)0733-9410(1991)117:9(1363).

[49] Saitoh M. On the performance of lumped parameter models with gyro-mass elements for the impedance function of a pile-group supporting a single-degree-of-freedom system. Earthq Eng Struct Dyn 2012;41:623-41.

[50] Saitoh M. Simple model of frequency-dependent impedance functions in soil-structure interaction using frequency-independent elements. J Eng Mech 2007;133:1101-14.

[51] Wang Jue, Ding Zhou WL, Wang S. Nested lumped-parameter model for foundation with strongly frequency-dependent impedance. J Earthq Eng 2016;20:975-91. http://dx.doi.org/10.1080/13632469.2015.1109568.

[52] CEA. Cast3m. 2001, http://www-cast3m.cea.fr/.

[53] Hardin B, Drnevich V. Shear modulus and damping in soils: Design equations and curves. J Soil Mech Found Div 1972;98. http://dx.doi.org/10.1061/JSFEAQ.0001760.

[54] Nohra C. Performance des Sols Argileux Renforcés par Inclusions Rigides sous Chargements Sismiques : Modélisation en Centrifugeuse et Numérique. (English title: Seismic Performance of Clay Soils Reinforced by Rigid Inclusions: Centrifuge and Numerical Modeling). (Ph.D. thesis), Université Gustave Eiffel; 2025, https://theses.hal.science/tel-05416397.

[55] IREX. Institute for Applied Research and Experimentation in Civil Engineering. ASIRI national project: Recommandations pour la conception, le calcul, l'exécution et le contrôle des ouvrages sur sols améliorés par inclusions rigides verticales. Presses des Ponts; 2012.

[56] Nohra C, Sorino-Camelo C, Escoffier S, Li Z, Thorel L. Experimental study of the kinematic interaction between a cohesive soil and a group of rigid inclusions through dynamic centrifuge tests. Jpn. Geotech. Soc. Spec. Publ. 2024;10(25):936-41. http://dx.doi.org/10.3208/jgssp.v10.OS-14-05.

[57] Coronado C, Gidwani N. Calculation of dynamic impedance of foundations using finite element procedures. San Francisco: Bechtel Power Corporation Nuclear security & Environmental; 2016.