# A fully-resolved micromechanical simulation of piping erosion during a suction bucket installation

Samuel Kemmler, Pablo Cuéllar, Antoni Artinov, Li-Hua Luu, Abbas Farhat, Pierre Philippe, Christoph Rettinger, Harald Köstler

► To cite this version:

Samuel Kemmler, Pablo Cuéllar, Antoni Artinov, Li-Hua Luu, Abbas Farhat, et al. A fully-resolved micromechanical simulation of piping erosion during a suction bucket installation. Computers and Geotechnics, 2025, 186, pp.107375. 10.1016/j.compgeo.2025.107375. hal-05228429

HAL Id: hal-05228429

https://hal.inrae.fr/hal-05228429v1

Submitted on 28 Aug 2025

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# A fully-resolved micromechanical simulation of piping erosion during a suction bucket installation

Samuel Kemmler<sup>a,b,D</sup>,*, Pablo Cuéllar<sup>c,a D</sup>, Antoni Artinov<sup>a D</sup>, Li-Hua Luu<sup>d</sup>, Abbas Farhat<sup>d</sup>, Pierre Philippe<sup>d D</sup>, Christoph Rettinger<sup>b D</sup>, Harald Köstler<sup>b D</sup>

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

## ARTICLE INFO

Dataset link: https://zenodo.org/records/14170253

MSC:
74L10

Keywords:

Micromechanical simulation
Fluid-solid coupling
Piping erosion
Suction bucket foundation
Offshore wind support structure
High-performance computing

## ABSTRACT

Granular fluidization phenomena such as piping erosion represent a challenge to the delicate installation process of offshore suction bucket foundations. A detailed analysis of the complex conditions in terms of soil composition, soil state, and foundation installation parameters that may lead to piping can be very demanding, if at all possible, solely by experimental means or using macroscopic continuum-based seabed models. The present paper presents a fully-resolved fluid-coupled micromechanical approach for a three-dimensional numerical simulation of the installation process of a suction bucket using the lattice Boltzmann method and discrete element method. The developed model is validated using well-established benchmarks and calibrated by means of experimental data from physical model tests on relevant scenarios focusing on the local fluidization of fixed embedded suction buckets as well as on the suction-driven installation of unrestrained buckets. The qualitative and quantitative agreement with the experimental data both endorse the proposed methodology and highlight the physical soundness of the obtained results. Thereby, the paper shows that three-dimensional analyses of relevant local scenarios at a real scale with little macromechanical model assumptions are feasible.

## 1. Introduction

Suction-aided intermediate foundations in their various forms (caissons or buckets, suction piles, suction anchors, etc.) have become increasingly popular in the past decades as a practical method to install offshore (sub)structures to the seabed (Kay et al., 2021; Sturm, 2017; Cotter, 2009). In particular, the multi-bucket jacket foundations show high potential for lightweight design of substructures for offshore wind turbines (Penner et al., 2020), see Fig. 1.

Such foundations are installed using the self-weight of the structure and induced hydraulic under-pressure due to suction inside the bucket, which combined generate a downward pressure over the bucket lid and effectively drives the bucket into the seabed. In the presence of cohesionless soils with significant permeability, like sandy soils or coarser grain size, the imposed pressure difference also creates a seepage flow through the soil into the bucket that alters the effective stress state in the soil matrix, thereby reducing the penetration resistance of the bucket. This installation approach is often faster, more cost-efficient,

reversible, and environmentally friendly, e.g., very low noise emissions, when compared to conventional methods like dynamic blow-driving of monopiles.

In general, the installation process is often the most delicate and dangerous phase in the bucket's service lifetime. Apart from the various structural threats derived from the high slenderness of the bucket wall, such as buckling phenomena or large plastic deformations upon boulder collision, there are several critical conditions that may arise during the suction-aided installation associated with the granular fluidization of the soil. From the so-called “rat-holing” at shallow embedment to localized events of piping erosion, the granular fluidization phenomena may lead to an effective breach of the bucket's hydraulic seal, thereby preventing the generation of differential pressures, eventually causing an altogether failure of the bucket installation process (Sturm, 2017). The broad range of conditions leading to granular fluidization during installation is complex and not sufficiently understood yet, whereby not only the soil conditions (heterogeneity, cementation, etc.) and state parameters (stress state and relative density) but also the practical

<footnote><sup>a</sup> Division 7.2 “Buildings and Structures”, Bundesanstalt für Materialforschung und -prüfung (BAM), Unter den Eichen 87, Berlin, 12205, Berlin, Germany</footnote>

<footnote><sup>b</sup> Chair for System Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, Erlangen, 91058, Bavaria, Germany</footnote>

<footnote><sup>c</sup> Referat 62 “Dams and Hydraulic constructions”, Bayerisches Landesamt für Umwelt, Bürgermeister-Ulrich-Straße 160, Augsburg, 86179, Bavaria, Germany</footnote>

<footnote><sup>d</sup> INRAE, Aix Marseille Univ., RECOVER, 3275 Route Cézanne, Aix-en-Provence, 13100, Provence-Alpes-Côte d'Azur, France</footnote>

<footnote>* Corresponding author at: Division 7.2 “Buildings and Structures”, Bundesanstalt für Materialforschung und -prüfung (BAM), Unter den Eichen 87, Berlin, 12205, Berlin, Germany.</footnote>

<footnote>E-mail address: samuel.kemmler@fau.de (S. Kemmler).</footnote>
Nomenclature

<table><tr><td>α</td><td>Calibration factor</td></tr><tr><td>q̄</td><td>Opposite lattice direction of q</td></tr><tr><td>δ<sub>ij,t</sub></td><td>Accumulated relative tangential motion between particles</td></tr><tr><td>c<sub>q</sub></td><td>Lattice velocity</td></tr><tr><td>f<sup>ext</sup></td><td>Force density on fluid</td></tr><tr><td>F<sub>p,i</sub></td><td>Total force on particle</td></tr><tr><td>F<sup>co</sup><sub>p,i</sub></td><td>Contact force on particle</td></tr><tr><td>F<sup>ext</sup><sub>p,i</sub></td><td>External force on particle</td></tr><tr><td>F<sup>hyd</sup><sub>p,i</sub></td><td>Hydrodynamic force on particle</td></tr><tr><td>F<sup>co</sup><sub>ij,n</sub></td><td>Normal part of particle contact force</td></tr><tr><td>F<sup>co,SD</sup><sub>ij,t</sub></td><td>Tangential part of spring-dashpot force</td></tr><tr><td>F<sup>co</sup><sub>ij,t</sub></td><td>Tangential part of particle contact force</td></tr><tr><td>n<sub>i,j</sub></td><td>Normal vector between particles</td></tr><tr><td>T<sub>p,i</sub></td><td>Total torque on particle</td></tr><tr><td>T<sup>co</sup><sub>p,i</sub></td><td>Contact torque on particle</td></tr><tr><td>T<sup>hyd</sup><sub>p,i</sub></td><td>Hydrodynamic torque</td></tr><tr><td>U<sub>f</sub></td><td>Fluid velocity</td></tr><tr><td>U<sub>p,i</sub>(x, t)</td><td>Velocity of particle i at position x</td></tr><tr><td>U<sub>p,i</sub></td><td>Particle velocity</td></tr><tr><td>U<sup>cp</sup><sub>ij,n</sub></td><td>Normal component of the relative velocity between particles</td></tr><tr><td>U<sup>cp</sup><sub>ij,t</sub></td><td>Tangential component of the relative velocity between particles</td></tr><tr><td>x</td><td>Lattice cell center</td></tr><tr><td>x<sub>i</sub></td><td>Cells intersecting with particle i</td></tr><tr><td>x<sub>p,i</sub></td><td>Particle position</td></tr><tr><td>x<sup>cp</sup><sub>i,j</sub></td><td>Contact point between particles</td></tr><tr><td>Δt</td><td>Time step size</td></tr><tr><td>Δx</td><td>Fluid cell size</td></tr><tr><td>Δx</td><td>Smagorinsky filter length scale</td></tr><tr><td>δ<sub>ij,n</sub></td><td>Penetration depth between particles</td></tr><tr><td>ε<sub>i</sub>(x, t)</td><td>Solid volume fraction of particle</td></tr><tr><td>Λ</td><td>Two-Relaxation-Time (TRT) magic number</td></tr><tr><td>λ<sub>+</sub>, λ<sub>-</sub></td><td>TRT relaxation parameters</td></tr><tr><td>μ</td><td>Dynamic viscosity</td></tr><tr><td>μ<sub>p</sub></td><td>Particle friction coefficient</td></tr><tr><td>ν</td><td>Kinematic viscosity</td></tr><tr><td>ν<sub>e</sub></td><td>Effective viscosity</td></tr><tr><td>ν<sub>p</sub></td><td>Poisson's ratio</td></tr><tr><td>ν<sub>t</sub></td><td>Eddy viscosity</td></tr><tr><td>ρ<sub>0</sub></td><td>Mean fluid density</td></tr><tr><td>ρ<sub>f</sub></td><td>Fluid density</td></tr><tr><td>ρ<sub>p</sub></td><td>Particle density</td></tr><tr><td>τ</td><td>Relaxation time</td></tr><tr><td>τ<sub>e</sub></td><td>Effective relaxation time coefficient</td></tr><tr><td>f̃<sub>q</sub></td><td>Post-collision Particle Distribution Functions (PDFs)</td></tr><tr><td>A</td><td>Surface area</td></tr><tr><td>B(x, t)</td><td>Total solid weighting factor</td></tr><tr><td>B<sub>i</sub>(x, t)</td><td>Solid weighting factor of particle i</td></tr><tr><td>C1 – C5</td><td>Suction bucket cases</td></tr><tr><td>C</td><td>Dimensionless normalized force</td></tr><tr><td>C<sub>q</sub></td><td>Collision operator</td></tr><tr><td>C<sub>S</sub></td><td>Smagorinsky constant</td></tr><tr><td>c<sub>S</sub></td><td>Lattice speed of sound</td></tr></table>

<table><tr><td>D</td><td>Diameter</td></tr><tr><td>d</td><td>Simulation domain depth</td></tr><tr><td>D<sub>p</sub></td><td>Particle diameter</td></tr><tr><td>E1 – E4</td><td>Settling sphere experiments</td></tr><tr><td>e</td><td>Bucket embedment into the soil</td></tr><tr><td>e<sub>dry</sub></td><td>Coefficient of restitution</td></tr><tr><td>F<sup>suction</sup></td><td>Suction force on bucket</td></tr><tr><td>F<sub>q</sub></td><td>Forcing operator on fluid</td></tr><tr><td>f<sub>q</sub></td><td>PDFs</td></tr><tr><td>f<sup>+</sup><sub>q</sub>, f<sup>-</sup><sub>q</sub></td><td>Symmetric and anti-symmetric parts of PDFs</td></tr><tr><td>f<sup>eq,+</sup><sub>q</sub>, f<sup>eq,-</sup><sub>q</sub></td><td>Symmetric and anti-symmetric parts of PDFs equilibrium</td></tr><tr><td>f<sup>eq</sup><sub>q</sub></td><td>PDFs equilibrium</td></tr><tr><td>g</td><td>Gravitational acceleration</td></tr><tr><td>h</td><td>Simulation domain height</td></tr><tr><td>i, j</td><td>Particle indices</td></tr><tr><td>k<sub>n</sub>, d<sub>n</sub></td><td>Normal stiffness and damping coefficients</td></tr><tr><td>k<sub>t</sub>, d<sub>t</sub></td><td>Tangential stiffness and damping coefficients</td></tr><tr><td>L</td><td>Normalized particle-wall distance</td></tr><tr><td>m<sub>p,i</sub></td><td>Particle mass</td></tr><tr><td>m<sub>ij,eff</sub></td><td>Particle effective mass</td></tr><tr><td>P</td><td>Pressure</td></tr><tr><td>Q</td><td>Suction rate</td></tr><tr><td>q</td><td>Lattice direction</td></tr><tr><td>R</td><td>Radius</td></tr><tr><td>Re</td><td>Reynolds number</td></tr><tr><td>Re<sub>s</sub></td><td>Shear Reynolds number</td></tr><tr><td>S<sub>ij</sub></td><td>Filtered strain rate tensor</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>T<sub>C</sub></td><td>Particle collision time</td></tr><tr><td>u<sub>p</sub></td><td>Particle settling velocity</td></tr><tr><td>u<sup>suction</sup></td><td>Suction velocity</td></tr><tr><td>u<sub>s</sub></td><td>Maximum particle settling velocity</td></tr><tr><td>ω</td><td>Simulation domain width</td></tr><tr><td>ω<sub>q</sub></td><td>Lattice weights</td></tr><tr><td>z</td><td>Particle coordination number</td></tr><tr><td>Ω<sub>p,i</sub></td><td>Rotational particle velocity</td></tr></table>

installation procedure, e.g., depth-dependent suction rates, may play a critical role.

In the past, physical model tests of bucket installation in enhanced-gravity conditions, e.g., in geotechnical centrifuges, have provided very useful data and a unique insight into the installation phenomena under different stress conditions (Ragni et al., 2020). However, such tests are generally very costly, difficult, and often at the limits of feasibility due to scaling effects, while the quantities of interest are also challenging to investigate, if at all measurable.

As a cost-efficient alternative, non-resolved coupled Computational Fluid Dynamics (CFD)-Discrete Element Method (DEM) simulations can provide some insights into the suction bucket installation (Yang et al., 2023; Zhang et al., 2024) and granular fluidization phenomena (Rettinger and Rüde, 2018). However, it is challenging to calibrate the non-resolved approaches to produce physically sound results. In contrast, the computationally demanding fully-resolved approaches allow very detailed and physically consistent analyses (Rettinger and Rüde, 2022) and have been shown capable of reproducing, both qualitatively and quantitatively, the complex physics of immersed granular
<figure><img src="image_4.png" /><figcaption>Fig. 1. Multi-bucket jacket foundation at the Borkum Riffgrund 2 installation (image gently provided by NGI/Ørsted).</figcaption></figure>

media with high accuracy (Schwarzmeier et al., 2023). As a con-
sequence, these approaches have often been applied to geotechnical
problems such as seepage failure and erosion (Froio et al., 2019; Tran
et al., 2017; Fukumoto et al., 2021; Fukumoto and Ohtsuka, 2018;
Benseghier et al., 2020; Kemmler et al., 2025a). Using efficient imple-
mentations, modern supercomputers can handle very large domains,
even for computationally demanding fully-resolved approaches (Kemm-
ler et al., 2025b). DEM simulations are a well-established technique
in geotechnical engineering to improve the understanding of complex
phenomena, such as pile installation by plain and rotary jacking in
soils (Cerfontaine et al., 2023) and cone penetration tests (Ciantia et al.,
2016).

The goal of the present work is to develop a comprehensive three-
dimensional numerical model for the installation of a suction bucket
foundation and compare it with experimental measurements consid-
ering fixed and unrestrained bucket conditions. The model aims to
simulate and analyze different suction conditions, in particular piping
erosion, to improve the understanding of the relationship between
the aforementioned broad range of conditions leading to granular
fluidization. Thereby, a fully-resolved micromechanical approach will
be applied to obtain detailed insights into the triggering conditions of
piping.

## 2. Numerical modeling

This section provides background information based on the works
published in Kemmler et al. (2025b) and Rettinger and Rüde (2017,
2022). Generally, fully-resolved coupled fluid-particle simulations con-
sist of three components: fluid dynamics, particle physics, and fluid-
particle coupling. In this context, “fully-resolved coupled fluid-particle”
refers to simulations with both a fluid and solid phase coupled through
momentum exchange between both phases. The fluid grid spacing is
significantly smaller than the particle diameter, which resolves flow
features around individual particles, allowing the accurate compu-
tation of hydrodynamic interaction for each particle, in contrast to
underresolved or continuum-based macroscopic approaches. Fig. 2
illustrates exemplary a fully-resolved coupled fluid-particle simulation.

The mathematical model applied in the present work is based on the
following main assumptions:

1. Purely-frictional granular soil: Particle-particle interactions are solely governed by unilateral frictional constraints, i.e., no inter-granular cohesion, adhesion, or cementation.

2. Perfectly spherical particles: Each particle of the soil is assumed to be a perfect sphere, neglecting shape effects such as angularity or elongation.

3. Linear contact elasticity: Particle contact forces are estimated with a Hookean elasticity model, in contrast to non-linear approaches like the Hertz-Mindlin contact model.

4. No particle deformations: All particles are assumed to be perfectly rigid, i.e., no elastic deformations nor plastic yielding, crushing, or breakage are considered. Only interpenetration between intact spherical particles in contact.

5. Incompressible fluid: Mach numbers are assumed to remain sufficiently low ($Ma < 0.1$) to minimize compressibility effects.

The mathematical model's assumptions certainly impact the results and impose limitations on its applicability to real-world scenarios where conditions deviate from these hypotheses. Regarding the solid phase, the adopted model is well suited for simulating purely frictional granular materials such as sands or gravels, whereas other types of soils, e.g., fine-grained soils, would require the introduction of additional intergranular forces for the consideration of cohesion, adhesion, cementation, etc., using suitable models (Potyondy and Cundall, 2004; Delenne et al., 2004; Benseghier et al., 2020). Besides, the assumption of spherical unbonded particles naturally precludes the appearance of specific intergranular effects such as interlocking or anisotropy due to particle shape, e.g., angularity, elongation, concavity, which are out of the scope of this paper. Concerning the constitutive behavior, all particles are considered perfectly rigid here, i.e., no changes in shape or size, while the visco-elastic rheological elements are just a convenient means to compute interaction forces. Therefore, particle deformation or degradation processes like attrition, breakage, crushing, etc., are also disregarded here, which appears, however, to be a reasonable assumption for the practical scenarios envisaged here, i.e., shallow soil layers with geostatic stress levels well below 100 kPa and little, if any, grain-crushing. Finally, the assumption of incompressible fluid hinders the simulation of compressible flow effects, such as shock waves, though these are not expected to be relevant in the current context.

This study focuses only on the short-term transient effects tak-
ing place during the installation of the suction foundation, which
are governed by the suction pressure applied within the bucket. The
long-term behavior and stability of the structure are, therefore, out
of scope here. Once the installation to target depth is achieved, no
further suction pressures are applied to the bucket anymore, and the
drained behavior becomes dominant in the long term. Nevertheless,
cyclic loading from wind and waves in multi-bucket foundations can
also induce over and under-pressures inside the bucket, so the present
model could also provide valuable insights for short-term scenarios
involving few load cycles. In any case, long-term cumulative effects
typically involve millions of load cycles, which are not tractable within
the modeling framework adopted here, considering the computation
cost. The model's explicit time-stepping scheme requires indeed a fine
temporal resolution, restricting simulations to short-term events on the
order of minutes. As a result, the methodology is only suitable for
capturing shorter effects, e.g., a few cycles for typical load frequencies
in the range of 0.1 Hz.

2.1. *Lattice Boltzmann method*

The Lattice Boltzmann Method (LBM) is employed for the fluid
dynamics simulation as an alternative to conventional Navier-Stokes
<figure><img src="image_5.png" /><figcaption>Fig. 2. Schematic representation of a two-dimensional fully-resolved coupled fluid-particle simulation with hydrodynamic force $F_{p,i}^{\text{hyd}}$ and torque $T_{p,i}^{\text{hyd}}$ acting on the particle i. Cells can either be fluid (outside particle, blue), solid (inside particle, orange), or partially saturated (intersecting with particle surface, light orange).</figcaption></figure>

solvers (Krüger et al., 2017). The LBM involves 27 PDFs $f_q$ with $q \in \{0, ..., 26\}$ for each cell of a three-dimensional Cartesian grid. Each $f_q$ is connected to a lattice velocity $c_q$. The cell is then updated based on the Boltzmann equation, split into the collision and streaming steps (Krüger et al., 2017). The cell-local collision step relaxes the PDFs towards their thermodynamic equilibrium state. Considering the lattice cell with center $\mathbf{x}$ at the time $t$, the collision step can be written as:

$$ \tilde{f}_q(\mathbf{x}, t) = f_q(\mathbf{x}, t) + C_q(\mathbf{x}, t) + F_q(\mathbf{x}, t), \quad (1) $$

with the collision operator $C_q$ and the forcing operator $F_q$. The streaming step propagates the post-collision PDFs $\tilde{f}_q$ to neighboring cells and is given as:

$$ f_q(\mathbf{x} + \Delta t c_q, t + \Delta t) = \tilde{f}_q(\mathbf{x}, t), \quad (2) $$

where $\Delta t$ is conveniently chosen as the time unit in the actual numerical implementation of LBM. The cell-local density $\rho_f$ and velocity $U_f$ are calculated using the moments of the PDFs:

$$ \rho_f(\mathbf{x}, t) = \sum_q f_q(\mathbf{x}, t), \quad (3) $$

$$ U_f(\mathbf{x}, t) = \frac{1}{\rho_0} \sum_q f_q(\mathbf{x}, t) c_q + \frac{\Delta t}{2\rho_0} f^{\text{ext}}, \quad (4) $$

using the mean fluid density $\rho_0 = 1$ and the constant external force density $f^{\text{ext}}$ (Ladd and Verberg, 2001). The TRT collision operator (Ginzburg and Verhaeghe, 2008) is given using two independent relaxation parameters $\lambda_+$, $\lambda_-$ in the range $(-2, 0)$:

$$ C_q^{\text{TRT}}(\mathbf{x}, t) = \lambda_+ \left[ f_q^+(\mathbf{x}, t) - f_q^{\text{eq,+}}(\rho_f, U_f) \right] + \lambda_- \left[ f_q^-(\mathbf{x}, t) - f_q^{\text{eq,-}}(\rho_f, U_f) \right], \quad (5) $$

with the PDFs and their equilibrium values being split into their symmetric and anti-symmetric parts:

$$ f_q^+ = \frac{1}{2}(f_q + f_{\bar{q}}), \quad f_q^- = \frac{1}{2}(f_q - f_{\bar{q}}), \quad (6) $$

with $\bar{q}$ being the inverse direction of $q$. $\lambda_+$ is related to the relaxation time $\tau$, i.e., $\lambda_+ = -\Delta t/\tau$. The relaxation time $\tau$ is linked to the kinematic viscosity $\nu$ by:

$$ \nu = (\tau - \Delta t/2)c_s^2, \quad (7) $$

with the lattice speed of sound $c_s = 1/\sqrt{3}$. $\Lambda$ relates the second relaxation rate $\lambda_-$ to the first one $\lambda_+$ through:

$$ \Lambda = \left( \frac{1}{2} + \frac{1}{\lambda_+} \right) \left( \frac{1}{2} + \frac{1}{\lambda_-} \right). \quad (8) $$

$\Lambda$ is commonly set to 3/16 (Khirevich et al., 2015). The equilibrium is defined as (Krüger et al., 2017):

$$ f_q^{\text{eq}}(\rho_f, U_f) = w_q \rho_f \left( 1 + \frac{c_q \cdot U_f}{c_s^2} + \frac{(c_q \cdot U_f)^2}{2c_s^4} - \frac{U_f \cdot U_f}{2c_s^2} \right), \quad (9) $$

with the lattice weights $w_q$ given as (Krüger et al., 2017):

<table><thead><tr><th>Lattice velocity $c_q^\mathrm{T}$</th><th>Lattice weight $w_q$</th></tr></thead><tbody><tr><td>(0, 0, 0)</td><td>8/27</td></tr><tr><td>(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)</td><td>2/27</td></tr><tr><td>(±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1)</td><td>1/54</td></tr><tr><td>(±1, ±1, ±1)</td><td>1/216</td></tr></tbody></table> \quad (10)

The forcing operator incorporates external forces and is defined as:

$$ F_q(\mathbf{x}, t) = \Delta t w_q \left( \frac{c_q - U_f}{c_s^2} + \frac{c_q \cdot U_f}{c_s^4} \cdot c_q \right) \cdot f^{\text{ext}}. \quad (11) $$

Furthermore, a subgrid-scale model is introduced to account for the effect of unresolved turbulence typically occurring due to high fluid velocities in the pores of the seabed. The (Smagorinsky, 1963) Large Eddy Simulation (LES) model modifies the computation of the first relaxation rate $\lambda_+$ by introducing the dimensionless effective relaxation time coefficient $\tau_e$ (Fukumoto et al., 2021):

$$ \lambda_+ = -\Delta t / \tau_e. \quad (12) $$

Using the effective viscosity $\nu_e = \nu + \nu_t$, the effective relaxation time coefficient $\tau_e$ can be obtained using Eq. (7). The eddy viscosity $\nu_t$ is calculated using the Smagorinsky constant $C_S$, the filter length scale $\Delta_x$, and the filtered strain rate tensor $S_{ij}$:

$$ \nu_t = (C_S \Delta_x)^2 |S_{ij}|. \quad (13) $$

In accordance with the literature, we use $C_S = 0.1$ (Yu et al., 2005) and $\Delta_x = 1$ (Fukumoto et al., 2021). The filtered strain rate tensor $S_{ij}$ is computed as:

$$ |S_{ij}| = - \frac{1}{2 \rho_f c_s^2 \tau_e} \sqrt{2 \sum_{i,j} Q_{ij} Q_{ij}}, \quad (14) $$

$$ Q_{ij} = \sum_q c_{q,i} c_{q,j} (f_q - f_q^{\text{eq}}). \quad (15) $$

## 2.2. Particle dynamics

The DEM (Cundall and Strack, 1979) is used to model the behavior of the particles. The total force $F_{p,i}$ and torque $T_{p,i}$ on a particle $i$ comprises of the following contributions:

$$ F_{p,i} = F_{p,i}^{\text{co}} + F_{p,i}^{\text{hyd}} + F_{p,i}^{\text{ext}}, \quad (16) $$

$$ T_{p,i} = T_{p,i}^{\text{co}} + T_{p,i}^{\text{hyd}}, \quad (17) $$

with the contact force $F_{p,i}^{\text{co}}$ and torque $T_{p,i}^{\text{co}}$, the hydrodynamic force $F_{p,i}^{\text{hyd}}$ and torque $T_{p,i}^{\text{hyd}}$, and the external force $F_{p,i}^{\text{ext}}$. The particle movements are obtained by integrating the equations of motion.

### 2.2.1. Particle contact handling

The contact interaction between particles $i$ and $j$ is modeled using a linear spring-dashpot model according to Cundall and Strack (1979). Thereby, the resulting contact force $F_{p,i}^{\text{co}}$ and torque $T_{p,i}^{\text{co}}$ acting on particle $i$ are given as:

$$ F_{p,i}^{\text{co}} = \sum_{j,j \neq i} (F_{ij,n}^{\text{co}} + F_{ij,t}^{\text{co}}), \quad (18) $$
$$T_{p,i}^{\text{co}} = \sum_{j,j \neq i} (\mathbf{x}_{ij}^{\text{cp}} - \mathbf{x}_{p,i}) \times F_{ij,t}^{\text{co}}, \qquad (19)$$

with the normal part of the contact force $F_{ij,n}^{\text{co}}$ at particle position $\mathbf{x}_{p,i}$ being computed as:

$$F_{ij,n}^{\text{co}} = -k_n \delta_{ij,n} \mathbf{n}_{ij} - d_n U_{ij,n}^{\text{cp}}, \qquad (20)$$

$k_n$ and $d_n$ are the normal stiffness and damping coefficients, respectively. $\mathbf{n}_{ij}$ is defined as the normal vector between particles i and j, $\delta_{ij,n}$ is the penetration depth during the contact, and $U_{ij,n}^{\text{cp}}$ denotes the normal part of the relative velocity between the particles at the contact point $\mathbf{x}_{ij}^{\text{cp}}$. The tangential part of the contact force is defined analogously as:

$$F_{ij,t}^{\text{co,SD}} = -k_t \delta_{ij,t} - d_t U_{ij,t}^{\text{cp}}, \qquad (21)$$

introducing the tangential stiffness and damping coefficients $k_t$ and $d_t$, respectively. $U_{ij,t}^{\text{cp}}$ denotes the tangential component of the relative velocity between the particles at the contact point. The accumulated relative tangential motion between two particles is defined as:

$$\delta_{ij,t} = \int_{t_i}^{t} U_{ij,t}^{\text{cp}}(t')dt', \qquad (22)$$

where $t_i$ is the instance in time of the impact according to Rettinger and Rüde (2022). The Coulomb friction model limits the tangential part of the collision force and incorporates slipping of the two particles by bounding the tangential force component from above using the coefficient of friction $\mu_p$, and the normal part of the collision force (Biegert et al., 2017):

$$F_{ij,t}^{\text{co}} = \min(|F_{ij,t}^{\text{co,SD}}|, |\mu_p F_{ij,n}^{\text{co}}|)I_{ij}, \qquad (23)$$

with

$$t_{ij} = F_{ij,t}^{\text{co,SD}} / |F_{ij,t}^{\text{co,SD}}|. \qquad (24)$$

The accumulated relative tangential motion $\delta_{ij,t}$ must not increase once the two particle surfaces are slipping and is correspondingly reset (Luding, 2008). The normal stiffness and damping coefficients $k_n$ and $d_n$ of the linear spring-dashpot model are commonly derived from the collision time $T_C$ and the coefficient of restitution $e_{dry}$ that denotes the ratio of the normal velocity after and before a collision:

$$k_n = \frac{m_{ij,\text{eff}}(\pi^2 + \ln^2 e_{\text{dry}})}{T_C^2}, \qquad (25)$$

$$d_n = -\frac{2m_{ij,\text{eff}} \ln e_{\text{dry}}}{T_C}, \qquad (26)$$

where the collision time $T_C$ quantifies the duration of a collision event in the absence of fluid. In accordance with Rettinger (2023), we choose:

$$T_C = 2.31 \frac{D_p}{c_s}. \qquad (27)$$

The effective mass is defined as:

$$m_{ij,\text{eff}} = \begin{cases} \frac{m_{p,i}m_{p,j}}{m_{p,j}+m_{p,i}}, & \text{sphere-sphere,} \\ m_{p,i}, & \text{sphere-wall.} \end{cases} \qquad (28)$$

The tangential stiffness and damping coefficients can be computed from the normal ones via the following relations (Rettinger, 2023):

$$k_t = \kappa_p k_n, \qquad (29)$$

$$d_t = \sqrt{\kappa_p} d_n, \qquad (30)$$

with

$$\kappa_p = \frac{2(1 - \nu_p)}{2 - \nu_p}, \qquad (31)$$

using the Poisson's ratio $\nu_p$.

### 2.2.2. Integration of the particle properties

The particle position $\mathbf{x}_{p,i}(t+\Delta t_p)$ and velocity $\mathbf{U}_{p,i}(t+\Delta t_p)$ are updated by solving the Newton–Euler equations of motion (Cundall and Strack, 1979) using the velocity Verlet integrator (Swope et al., 1982).

### 2.3. Fully resolved fluid–particle coupling method

The task of the coupling is to exchange momentum between the fluid and the solid phases. The Partially Saturated Cells Method (PSM) is used for a fully resolved fluid–particle coupling. The idea behind the PSM is to modify the LBM collision step from Eq. (1) as follows:

$$\tilde{f}_q(\mathbf{x}, t) = f_q(\mathbf{x}, t) + C_q^{\text{PSM}}(\mathbf{x}, t) + (1 - B(\mathbf{x}, t))F_q(\mathbf{x}, t), \qquad (32)$$

where $B(\mathbf{x}, t)$ is the total solid weighting factor, which is the sum over the individual solid weighting factors $B_i(\mathbf{x}, t)$ of all particles i overlapping with the fluid cell $\mathbf{x}$, i.e., $B(\mathbf{x}, t) = \min(\sum_i B_i(\mathbf{x}, t), 1)$. $B_i(\mathbf{x}, t)$ is directly related to the solid volume fraction $\epsilon_i(\mathbf{x}, t)$, i.e., the fraction of a fluid cell $\mathbf{x}$ being partly or completely covered by particle i. In the present work, the solid weighting factor is (Noble and Torczynski, 1998):

$$B_i(\mathbf{x}, t) = \frac{\epsilon_i(\mathbf{x}, t)\left(\frac{\tau}{\Delta t} - \frac{1}{2}\right)}{(1 - \epsilon_i(\mathbf{x}, t)) + \left(\frac{\tau}{\Delta t} - \frac{1}{2}\right)}. \qquad (33)$$

The modified PSM collision operator $C_q^{\text{PSM}}$ used in Eq. (32) is given as (Wang et al., 2018):

$$C_q^{\text{PSM}}(\mathbf{x}, t) = (1 - B(\mathbf{x}, t))C_q^{\text{TRT}}(\mathbf{x}, t) + \sum_i B_i(\mathbf{x}, t)C_{q,i}^{\text{solid}}(\mathbf{x}, t). \qquad (34)$$

The solid collision operator $C_{q,i}^{\text{solid}}(\mathbf{x}, t)$ is only active for fluid cells with overlapping particles. In the present work, the solid collision operator is computed as (Noble and Torczynski, 1998):

$$C_{q,i}^{\text{solid}}(\mathbf{x}, t) = [f_q^{\text{eq}}(\rho_f, U_{p,i}(\mathbf{x}, t)) - f_q(\mathbf{x}, t)] + \left(1 - \frac{\Delta t}{\tau}\right)[f_q(\mathbf{x}, t) - f_q^{\text{eq}}(\rho_f, U_f)], \qquad (35)$$

where $f_q$ are the inverse PDFs of $f_q$. $U_{p,i}(\mathbf{x}, t)$ is the particle velocity with respect to the cell center $\mathbf{x}$:

$$U_{p,i}(\mathbf{x}, t) = U_{p,i}(t) + \Omega_{p,i}(t) \times [\mathbf{x} - \mathbf{x}_{p,i}(t)], \qquad (36)$$

using the particle's translational and rotational velocity $\mathbf{U}_{p,i}(t)$ and $\Omega_{p,i}(t)$, and the particle center of mass $\mathbf{x}_{p,i}(t)$. The hydrodynamic force $\mathbf{F}_{p,i}^{\text{hyd}}(t)$ and resulting torque $\mathbf{T}_{p,i}^{\text{hyd}}(t)$ acting from the fluid on particle i are computed as follows:

$$\mathbf{F}_{p,i}^{\text{hyd}}(t) = \frac{(\Delta x)^3}{\Delta t} \sum_{\mathbf{x}_i} [B_i(\mathbf{x}_i, t) \sum_q (C_{q,i}^{\text{solid}}(\mathbf{x}_i, t) c_q)], \qquad (37)$$

$$\mathbf{T}_{p,i}^{\text{hyd}}(t) = \frac{(\Delta x)^3}{\Delta t} \sum_{\mathbf{x}_i} [B_i(\mathbf{x}_i, t)(\mathbf{x}_i - \mathbf{x}_{p,i}) \times \sum_q (C_{q,i}^{\text{solid}}(\mathbf{x}_i, t)c_q)], \qquad (38)$$

where $\mathbf{x}_i$ are the centers of all cells intersecting with the particle.

### 2.4. Particle mapping

Jones and Williams (2017) addressed the problem that for the computation of $\epsilon_i(\mathbf{x}, t)$ for spheres, in general, no unique analytical solution exists. In their work, a linear approximation approach is introduced based on the analytical solution for a specific cell orientation relative to the sphere surface. They show that the linear approximation yields accurate results for varying particle orientations and is more computationally efficient than competing techniques like sub-division sampling. Therefore, it is implemented in the proposed model.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Schematic of the three tests: (left) sphere in Stokes flow; (middle) sphere in shear flow; (right) settling sphere.</figcaption></figure>

## 2.5. Notes on computational effort

Fully-resolved coupled fluid-particle simulations are computationally very intense and require efficient implementation and parallelization. The model proposed in this work is implemented in the massively parallel multiphysics framework WALBERLA (Bauer et al., 2021) (https://www.walberla.net/). The implementation and parallelization are introduced in Kemmler et al. (2025b) and involves a hybrid approach, i.e., the fluid simulation runs on Graphics Processing Units (GPUs), and the particle simulation on Central Processing Units (CPUs). A block-structured domain partitioning approach is used to achieve parallelism, i.e., the simulation domain is divided into multiple uniform blocks, which are exclusively assigned to one GPU each. The parallelization requires the communication of fluid and particle information between neighboring blocks. It was shown that the fluid simulation using the PSM, which is the computationally most demanding part of this fully-resolved approach, reaches high computational efficiency in terms of utilizing the main memory bandwidth (Kemmler et al., 2025b). Furthermore, due to pure nearest neighbor communication, WALBERLA has demonstrated good scalability up to thousands of GPUs (Holzer et al., 2021). Therefore, the model can run large-scale simulations up to real-size problems on thousands of GPUs on modern supercomputers. The code, simulation data, and corresponding scripts are available on Zenodo, see the data availability statement.

# 3. Model validation

This section provides a validation of the methodology by employing the following common simulation tests to validate the fluid-particle coupling for a single sphere illustrated in Fig. 3: a fixed sphere surrounded by Stokes flow, a fixed sphere exposed to shear flow, and a freely settling sphere under gravity. The tests and corresponding sections are based on a previous set of validation tests introduced in Rettinger and Rüde (2022).

## 3.1. Sphere in Stokes flow

The configuration considered here consists of a steady-state Stokes flow around an infinite array of static spheres. The simulation results are compared with a quasi-analytical solution. The domain is periodic in all directions and has a size of $2D_p \times 2D_p \times 2D_p$, with the particle diameter $D_p = 20$ fluid cells. The spherical particle is located in the center of the domain. The flow is driven by a constant external body force on the fluid, leading to an acceleration in the flow direction of $(1 \cdot 10^{-8}, 0, 0)$ lu/tu² in lattice units, where “lu” denotes the lattice cell size and “tu” the time step size. The reference value for the dimensionless normalized force is given as $C^{\text{ref}} = 2.842$ for this configuration (Sangani and Acrivos, 1982). The relative error between the simulated force and the reference value for different viscosities is shown in Fig. 4.

<figure><img src="image_7.png" /><figcaption>Fig. 4. Relative error of the drag force vs. fluid viscosity. Note that the viscosity is given in lattice units.</figcaption></figure>

The relative error starts at 3.9% for a kinematic viscosity of $0.0167 \text{ lu}^2/\text{tu}$ and increases to around 15.5% for a viscosity of $0.5 \text{ lu}^2/\text{tu}$. One can conclude that the accuracy of the numerically obtained relative error strongly depends on the viscosity of the fluid (smaller errors for lower viscosities). This observation aligns with the findings of Rettinger and Rüde (2017), who conclude that the PSM exhibits a viscosity-dependent behavior, and they report similar relative errors for this benchmark. For the suction bucket model presented below, a viscosity of $9.42 \cdot 10^{-3} \text{ lu}^2/\text{tu}$ is used, which corresponds to an error of 3.5%.

## 3.2. Sphere in shear flow

In this simulation test, the computational domain comprises a fixed sphere next to a resting plane in a uniform shear flow setup. This configuration enables the test of the model's capability to predict drag and lift forces on a sphere for different Reynolds numbers. As a reference, empirically obtained drag and lift coefficients taken from Zeng et al. (2009) are used. The domain is periodic in x- and y-direction and has a size of $48D_p \times 16D_p \times 8D_p$ with a sphere diameter of $D_p = 20$ fluid cells. The sphere is horizontally centered, i.e., located at position $(24D_p, 8D_p, LD_p)$ with a given dimensionless distance L from the bottom wall. In contrast to the bottom wall, which is at rest, the top wall is moving with a constant horizontal velocity $(0.1, 0, 0)$. The flow characteristic is described by the shear Reynolds number $Re_s$ as shown in Rettinger and Rüde (2022). The comparison between the drag and lift coefficients of the simulations and the reference solution Zeng et al. (2009) for different particle-wall distances and flow configurations is given in Table 1.

The simulation results show a minimal and maximal error of 2.0% and 15.6%, respectively. The error for the lift coefficient is higher
<table><caption>**Table 1**
Comparison of the drag and lift coefficients between the reference solution obtained in Zeng et al. (2009) and the present simulation results for different particle-wall distances L and flow configurations $Re_s$.</caption><thead><tr><th>L</th><th>Re<sub>s</sub></th><th>Coeff.</th><th>Ref.</th><th>Sim.</th><th>Rel. err. (%)</th></tr></thead><tbody><tr><td>0.505</td><td>1</td><td>Drag</td><td>44.88</td><td>41.71</td><td>7.1</td></tr><tr><td>0.505</td><td>1</td><td>Lift</td><td>3.54</td><td>3.10</td><td>12.4</td></tr><tr><td>0.505</td><td>25</td><td>Drag</td><td>3.53</td><td>3.60</td><td>2.0</td></tr><tr><td>0.505</td><td>25</td><td>Lift</td><td>0.88</td><td>0.83</td><td>5.7</td></tr><tr><td>1</td><td>1</td><td>Drag</td><td>35.77</td><td>32.31</td><td>9.7</td></tr><tr><td>1</td><td>1</td><td>Lift</td><td>2.32</td><td>2.06</td><td>11.2</td></tr><tr><td>1</td><td>25</td><td>Drag</td><td>2.81</td><td>2.70</td><td>3.9</td></tr><tr><td>1</td><td>25</td><td>Lift</td><td>0.45</td><td>0.38</td><td>15.6</td></tr></tbody></table>

<table><caption>**Table 2**
Fluid density $ρ_f$, dynamic viscosity $μ$, maximum particle settling velocity $u_t$, and resulting particle Reynolds number $Re$ for the four experiments E1-E4 from Ten Cate et al. (2002).</caption><thead><tr><th>Exp.</th><th>ρ<sub>f</sub> (kg/m<sup>3</sup>)</th><th>μ (Pa s)</th><th>u<sub>s</sub> (m/s)</th><th>Re</th></tr></thead><tbody><tr><td>E1</td><td>970</td><td>0.373</td><td>0.036</td><td>1.4</td></tr><tr><td>E2</td><td>965</td><td>0.212</td><td>0.057</td><td>3.9</td></tr><tr><td>E3</td><td>962</td><td>0.113</td><td>0.087</td><td>11.1</td></tr><tr><td>E4</td><td>960</td><td>0.058</td><td>0.122</td><td>30.3</td></tr></tbody></table>

than for the drag coefficient for all particle-wall distances and flow configurations. The error for the low $Re_s$ flow configurations is mostly larger than for the high $Re_s$ simulations.

### 3.3. Settling sphere

A settling sphere configuration is selected, replicating the fall of a single sphere under the influence of a gravity field within a surrounding fluid. Therefore, this comparison validates both the coupling and the time integration of the particle location and velocity. The temporal evolution of the sphere distance from the bottom wall and the particle settling velocity are compared with experimental data from Ten Cate et al. (2002) for four different silicon oils with varying viscosities. The experimental domain has the dimensions 0.1 m × 0.1 m × 0.16 m and contains a spherical particle with a diameter $D_p = 0.015$ m and a density of 1120 kg/m³. The initial distance between the sphere surface and the bottom wall is 8.25$D_p$. The characteristics of the four distinct experimental cases (E1-E4) with varying fluid properties are summarized in Table 2 together with the maximum sphere settling velocity measured in the experiments. In the simulation, the domain is discretized using 135 × 135 × 216 cells, resulting in a particle resolution of $D_p/Δx = 20.25$ cells. Fig. 5 shows the simulation and experimental data for the gap width between the particle and the bottom wall (normalized by the particle diameter) and the particle settling velocity vs. time. There is an excellent agreement between the simulation and the experimental results in all four cases regarding both sphere-wall distance and settling velocity. For higher viscosities, the deviation between numerical and experimental data is slightly higher than for lower viscosities. The relative errors of the maximum settling velocity for E1-E4 are 3.2%, 4.3%, 2.0%, and 0.2%, respectively. Hence, the coupling and time integration schemes accurately simulate the settling behavior at moderate Reynolds numbers.

## 4. Suction bucket installation process

The numerical study carried out in the present work is specifically devoted to the application case of suction bucket foundations during the installation of offshore wind turbines. The installation process is challenging due to a high risk of suction-induced destabilization of the foundation soil by piping erosion and appears particularly well suited to the practical testing of the coupled micro-mechanical model presented in Section 2. In this section, the numerical model of the suction bucket setup is described. An experimental study of suction-assisted bucket

settling in granular soil carried out specifically to produce benchmark data is presented. A comparison between numerical and experimental results and an in-depth analysis of the micromechanics are provided.

### 4.1. Numerical setup

The computational domain covers a representative sub-region of a suction bucket foundation, see Fig. 6. The granular soil is modeled using spherical particles, whereas a box represents the bucket wall segment. The domain is fully filled with fluid. The outflow is implemented using a velocity boundary condition mimicking the suction applied inside the bucket during the installation process. The outflow boundary condition is acting on the entire top right-hand side, in contrast to the experimental setup, which employs the suction via a pipe in the center of the bucket. However, applying the outflow on the entire top right-hand side has no significant impact on the system behavior: If the outflow boundary condition is applied to only a third of the top right-hand side with the velocity being three times as high, i.e., the suction mass is conserved, the differential pressure in a quasi-static case only differs by around 0.2%. The inflow is induced by a constant density/pressure boundary condition on the top left-hand side, reflecting hydro-static conditions outside the bucket. The domain is periodic in the out-of-plane direction in Fig. 6, i.e., the fluid and particles leaving the domain on one side re-enter on the opposite side. No-slip boundary conditions are applied on the remaining boundaries (left, top middle, right, and bottom), which are far enough away from the preferential flow path to have no significant influence on the system behavior (Saint-Venant's principle): Switching from no-slip to free-slip boundary conditions on the left and right walls in a quasi-static case changes the converged differential pressure only by around 1.9%. The bucket wall segment is handled by the PSM, similar to the spherical particles of the granular soil. The model, in principle, offers the possibility to activate lubrication correction terms. However, the present study focuses on the analysis of settled dense granular assemblies with most contacts in static conditions and relatively few particles approaching each other in the limit of vanishing interparticle gaps. Therefore, lubrication forces are assumed to be negligible, which was confirmed by comparing two setups, one with activated lubrication forces and one without, which resulted in no significant difference in the differential pressure. Thus, the lubrication correction is disabled in the present study for the sake of simplicity. The interested reader may find further details on the flow resolution in the limit of vanishing interparticle gaps in Rettinger and Rüde (2022). The simulation parameters specified in the remainder of this section are chosen to match the experimental setup described in Section 4.2. The lateral walls are $d = 2.1$ cm apart. The other two dimensions of the domain are $h \approx 13$ cm and $w \approx 11$ cm. The bucket wall segment can either be fixed or move unrestrained in the vertical direction due to a given force and is initially embedded by $e \approx 2$ cm into the granular soil. The differential pressure between the inside and outside of the bucket due to the applied suction velocity can be obtained from the density difference, which is monitored next to the bucket wall at 95% of the domain height, see Fig. 6. From the density difference, the differential pressure is computed using the following relation: $ΔP = c_s^2 \cdot Δρ_f$. The particle bed has a porosity of about 36% and consists of 10257 particles with a diameter of 3 mm (±0.3 mm). Since the experimental values are unknown, based on previous work, a friction coefficient of 0.16, a restitution coefficient of 0.926, and a Poisson's ratio of $v_p = 0.245$ are used (Tang et al., 2019). Note that LBM simulations are typically not performed in SI units but in lattice units. These lattice units are an artificial set of units and, therefore, require converting from SI units into lattice units and vice versa (Krüger et al., 2017). The simulation parameters in both lattice and SI units are given in Table 3. Note that the conversion factors for space $x$, time $t$, and density $ρ_f$ are explicitly given in the table from which the conversion factors for all other quantities can be derived. While the reference scales for $x$, $t$, and $ρ_f$ in SI units are given as m, s, and kg/m³, these are
<figure><img src="image_8.png" /><figcaption>Fig. 5. For E1-E4 and the corresponding simulation results: (left) normalized gap between the sphere and the bottom wall vs. time; (right) velocity of the sphere vs. time.</figcaption></figure>

<table><caption>**Table 3**
Simulation parameters in lattice and SI units.</caption><thead><tr><th>Parameter</th><th>Lattice units</th><th>SI units</th></tr></thead><tbody><tr><td>Domain size <i>w</i> × <i>d</i> × <i>h</i></td><td>740 lu × 140 lu × 870 lu</td><td>0.111 m × 0.021 m × 0.1305 m</td></tr><tr><td>Cell size Δ<i>x</i></td><td>1 lu</td><td>1.5 × 10<sup>-4</sup> m</td></tr><tr><td>Bucket slice size</td><td>67.34 lu × 140 lu × 336.69 lu</td><td>0.010 m × 0.021 m × 0.051 m</td></tr><tr><td>Particle diameter <i>D</i><sub>p</sub></td><td>18 lu to 22 lu</td><td>2.7 × 10<sup>-3</sup> m to 3.3 × 10<sup>-3</sup> m</td></tr><tr><td>Particle Density <i>ρ</i><sub>p</sub></td><td>2.25 <i>pu</i></td><td>2650 kg/m<sup>3</sup></td></tr><tr><td>Gravitational acceleration <i>g</i></td><td>6.54 × 10<sup>-6</sup> lu/tu<sup>2</sup></td><td>9.81 m/s<sup>2</sup></td></tr><tr><td>Fluid density <i>ρ</i><sub>f</sub></td><td>1 <i>pu</i></td><td>1179 kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity (kinematic) <i>ν</i></td><td>9.42 × 10<sup>-5</sup> lu<sup>2</sup>/tu</td><td>2.12 × 10<sup>-5</sup> m<sup>2</sup>/s</td></tr><tr><td>Runtime</td><td>1.5 × 10<sup>6</sup> tu</td><td>15 s</td></tr><tr><td>Time step size Δ<i>t</i></td><td>1 <i>tu</i></td><td>1 × 10<sup>-5</sup> s</td></tr><tr><td>Restitution coefficient <i>e</i><sub>dry</sub></td><td>0.926</td><td>0.926</td></tr><tr><td>Friction coefficient</td><td>0.16</td><td>0.16</td></tr><tr><td>Collision time <i>T</i><sub>C</sub></td><td>88 <i>tu</i></td><td>88 × 10<sup>-5</sup> s</td></tr><tr><td>Poisson's ratio <i>ν</i><sub>p</sub></td><td>0.245</td><td>0.245</td></tr><tr><td>Normal spring stiffness <i>k</i><sub>n</sub> (particle-particle)</td><td>7.99 <i>pu</i> lu<sup>3</sup> tu<sup>-2</sup></td><td>318 N m<sup>-1</sup></td></tr><tr><td>Normal spring stiffness <i>k</i><sub>n</sub> (particle-wall)</td><td>15.97 <i>pu</i> lu<sup>3</sup> tu<sup>-2</sup></td><td>636 N m<sup>-1</sup></td></tr><tr><td>Tangential spring stiffness <i>k</i><sub>t</sub> (particle-particle)</td><td>6.87 <i>pu</i> lu<sup>3</sup> tu<sup>-2</sup></td><td>273 N m<sup>-1</sup></td></tr><tr><td>Tangential spring stiffness <i>k</i><sub>t</sub> (particle-wall)</td><td>13.74 <i>pu</i> lu<sup>3</sup> tu<sup>-2</sup></td><td>547 N m<sup>-1</sup></td></tr><tr><td>Normal spring damping <i>d</i><sub>n</sub> (particle-particle)</td><td>10.94 <i>pu</i> lu<sup>3</sup> tu<sup>-1</sup></td><td>4.35 × 10<sup>-3</sup> N s m<sup>-1</sup></td></tr><tr><td>Normal spring damping <i>d</i><sub>n</sub> (particle-wall)</td><td>21.89 <i>pu</i> lu<sup>3</sup> tu<sup>-1</sup></td><td>8.71 × 10<sup>-3</sup> N s m<sup>-1</sup></td></tr><tr><td>Tangential spring damping <i>d</i><sub>t</sub> (particle-particle)</td><td>10.15 <i>pu</i> lu<sup>3</sup> tu<sup>-1</sup></td><td>4.04 × 10<sup>-3</sup> N s m<sup>-1</sup></td></tr><tr><td>Tangential spring damping <i>d</i><sub>t</sub> (particle-wall)</td><td>20.30 <i>pu</i> lu<sup>3</sup> tu<sup>-1</sup></td><td>8.08 × 10<sup>-3</sup> N s m<sup>-1</sup></td></tr></tbody></table>

denoted as lu, tu, and pu in lattice units, respectively. The remainder of this section summarizes the creation of the initial seabed conditions. First, a randomly packed particle bed is created by settling spheres onto the ground under the influence of a gravity field. No fluid is used yet. Second, particles that would overlap with the 2 cm embedded bucket wall segment are removed. Third, the bucket wall segment is added. Finally, the particles are settled toward the bucket wall segment. Following these pre-steps, fluid is added, and the coupled simulation can start.

## 4.2. Experimental setup

In the following, the experiments are presented, i.e., the material and sample preparation, the experimental setup, and the experimental tests selected to produce benchmark data. These experiments were carried out at the geomechanics laboratory of the RECOVER unit (INRAE, Aix-Marseille University) in Aix-en-Provence (France).

### 4.2.1. Material and sample preparation

The granular material used in the experiments consists of 3 mm mono-sized silicate glass beads. Suction flow is generated by a gear

pump with a maximum flow rate of 32 L/min. A more viscous liquid than water was used to reach higher pressure levels, making it easier for the bucket to penetrate the granular soil. The chosen liquid is a mixture of 74% glycerol and 26% distilled water, with a dynamic viscosity of around 0.024 Pa · s (±0.001 Pa · s), measured by a falling ball viscometer. A glass bead mass of 15 kg, supplemented by liquid to the extent of 30% of the total mass, was used. After full saturation of the granular phase, the sample is prepared by randomly stirring and folding in the glass beads to get rid of any memory effect on the granular phase from previous experiments. Once the whole granular sample has been prepared, a ruler is used to scrape off the excess beads, obtaining a smooth, flat surface.

### 4.2.2. Experimental setup

Figs. 7 and 8 show the experimental setup designed to reproduce the field installation of a suction bucket foundation at a small scale. This configuration enables the study of the fundamental behavior of bucket embedment, along with determining the critical flow rate and hydraulic pressure gradient required for successful installation at this reduced scale and with these model materials. The device consists of a cylindrical bucket with an internal diameter of *D* = 10 cm, see Fig. 8.
<figure><img src="image_9.png" /><figcaption>Fig. 6. (left) Schematic of an offshore wind turbine with a suction bucket foundation during installation; (right) three-dimensional numerical setup with grey solid particles, the dark bucket slice segment, and a 2D slice of the fluid velocity.</figcaption></figure>

The bucket can either be moved unrestrained along the vertical axis or fixed at a specific depth. The bucket is connected to a suction system powered by a gear pump. The bucket is initially embedded by $e = 2$ cm into the granular soil. The sample container and bucket are enclosed in a vessel, 57.5 cm by 52.7 cm in cross-section, partially filled with liquid. There are two different hydraulic circuits from the pump. The first one to the saturation inlet, highlighted in green in Fig. 8, saturates the sample at a small flow rate *before* starting the experiment. The second circuit, highlighted in red in Fig. 8, imposes a suction from the top outlet of the bucket and returns the liquid *via* a distributor back to the vessel. Furthermore, a counterweight has been added with a mass and pulley system to compensate, partially or totally, for the load applied by the bucket, valves, and pipe assembly. A camera is mounted in front of the vessel, recording the embedment of the bucket over time.

### 4.2.3. Experimental tests selected for comparison

The experiments carried out for the cross-comparison with numerical simulations comprise three cases where the bucket is fixed and two cases where the bucket is permitted for unrestrained movement along the vertical axis into the granular soil due to suction. The experimental conditions are given in Table 4. Fig. 9 presents selected frames of the observed behavior when increasing the suction rate in the granular sample with a fixed bucket for C1-C3. For the lowest flow rates of 1.4 L/min and 2.8 L/min (C1 and C2), no significant changes are detected at any of the three selected times, namely $t = 0$ s, 10 s, and 15 s. For the highest flow rate of 4.2 L/min (C3), a small heave starts at $t = 7$ s, eventually leading to a localized piping at $t = 8$ s, which can be seen in the frames at $t = 10$ s and $t = 15$ s. In Fig. 10, selected frames at $t = 0$ s, 10 s and 15 s are presented for C4 and C5, i.e., unrestrained bucket conditions. At the lowest flow rate of 1.6 L/min (C4), only very limited, yet detectable, movement of the bucket is found. In comparison, at the higher flow rate of 3.1 L/min (C5), the bucket displays much higher dynamics: Significant heave is visible in the granular sample, accompanied by a marked embedment of the bucket into the soil. The bucket embedment and pressure evolution measurements that were carried out during the five cases will be presented in Section 4.3.

<figure><img src="image_10.png" /><figcaption>Fig. 7. Experimental suction bucket embedment setup.</figcaption></figure>

### 4.3. Cross-comparison of experimental and numerical results

In this section, numerical results are compared with experimental data both for a fixed and unrestrained bucket.
<figure><img src="image_11.png" /><figcaption>Fig. 8. (left) Schematic of the experimental suction bucket embedment setup; (right) close-up image of the bucket displaying the position of the pressure ports.</figcaption></figure>

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

Suction rates *Q* and bucket freedom condition, i.e., fixed or unrestrained, of the experiments selected for cross-comparison.

<table><thead><tr><td>Case</td><td>Q (L/min)</td><td>Bucket</td></tr></thead><tbody><tr><td>C1</td><td>1.4</td><td>fixed</td></tr><tr><td>C2</td><td>2.8</td><td>fixed</td></tr><tr><td>C3</td><td>4.2</td><td>fixed</td></tr><tr><td>C4</td><td>1.6</td><td>unrestrained</td></tr><tr><td>C5</td><td>3.1</td><td>unrestrained</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Selected frames at *t* = 0 s, 10 s and 15 s of the fixed bucket over the granular sample for C1, C2 and C3 with increasing suction rates of 1.4 L/min, 2.8 L/min, and 4.2 L/min, respectively. The bed surface is indicated using dashed white lines.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Selected frames at *t* = 0 s, 10 s, and 15 s of the unrestrained bucket over the granular sample for C4 and C5 with increasing suction rates of 1.6 L/min and 3.1 L/min, respectively. The bed surface is indicated using dashed white lines.</figcaption></figure>

### 4.3.1. Determination of suction velocity and suction force

The experimental data provides the applied suction rate in liters per minute. In contrast, in the numerical simulation, the outflow boundary condition requires a velocity, i.e., meters per second. To convert the experimental suction rate *Q* (L/min) into the numerical suction velocity *u*<sub>suction</sub> (m/s), the following relationship is proposed:

$$u_{\text{suction}} = \frac{1}{1000 \text{ L/m}^3 \cdot 60 \text{ s min}^{-1}} \cdot \frac{Q}{A_{\text{bucket}}} \cdot \alpha, \quad (39)$$

with the experimental bucket cross-section area $A_{bucket} = \pi R_{bucket}^2$, the bucket radius $R_{bucket} = 5 \times 10^{-2}$ m, and a calibration factor $\alpha$, which has to be determined. This calibration factor accounts for the fact that the computational domain is a rectangular container and not a cylindrical bucket slice. A rectangular container and a cylindrical slice differ in shape and proportion as illustrated in Fig. 11. The wall segment depth *d* in the numerical setup is similar to the cylindrical slice of the experimental bucket. In contrast, the bed surface areas
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Zenithal perspective on (left) a cylindrical bucket slice; (right) the rectangular computational domain. *d*, *R<sub>bucket</sub>*, *A<sub>exp</sub>* and *A<sub>sim</sub>* denote the simulation domain depth, bucket radius, bed surface area of the cylindrical slice, and bed surface area of the simulation domain, respectively.</figcaption></figure>

<table><caption>**Table 5**

Calibration factor $\alpha$ for three experimental flow quantities to be preserved in the simulation.</caption><thead><tr><th>Exp. quantity preserved in sim.</th><th>α</th></tr></thead><tbody><tr><td>Flow rate underneath the bucket wall</td><td>0.5</td></tr><tr><td>Differential pressure for quasi-static case (empirically obtained)</td><td>0.63</td></tr><tr><td>Flow rate per bed surface area</td><td>1.0</td></tr></tbody></table>

$A_{exp}$ and $A_{sim}$ differ by a factor of around 2. Due to this geometric difference, not all flow quantities can be preserved at the same time. In Table 5, three possible values for the calibration factor $\alpha$, together with the corresponding preserved flow quantity, are provided. For $\alpha = 1.0$, the experimental suction rate is divided by the bucket surface area, resulting in the average experimental outflow velocity. If this velocity is applied as the outflow boundary condition in the numerical setup, the flow rate per bed surface area is the same in the experiment and simulation. However, since $A_{sim} \approx 2 \cdot A_{exp}$, the flow rate underneath the bucket wall segment is around twice as high in the simulation as in the experiment. For $\alpha = 0.5$, the average outflow velocity in the numerical setup is halved compared to the experimental setup. But since $A_{sim}$ is around two times $A_{exp}$ with a similar wall segment depth, the flow rate underneath the wall segment is preserved. Another important flow quantity is the differential pressure. A quasi-static case, i.e., small suction to avoid grain displacement, is used to determine $\alpha$ such that the differential pressure matches between the simulation and experiment. This can be achieved with the C1 experiment. This empirical approach yields $\alpha = 0.63$. Since an *a priori* estimate of how to choose $\alpha$ is not straightforward, the influence of $\alpha$ on the numerical results is investigated for different suction rates. The calibration factor $\alpha$ effectively serves as a scaling parameter for the numerical suction velocity. It addresses the geometrical differences between the simulation domain and the experimental setup, which prevent a direct derivation of numerical suction velocity from the experimental suction rate. Three values for $\alpha$ are proposed: 0.5, 0.63, and 1.0. The values 0.5 and 1.0 arise from theoretical considerations based on geometrical differences. In contrast, 0.63 is empirically obtained by calibrating the numerical suction velocity in a quasi-static flow regime to match the differential pressure observed in a corresponding experiment. In Section 4.3.2, the three proposed values are assessed by comparing numerical and experimental results to determine the most suitable calibration factor used by subsequent analyses. Since $\alpha$ directly scales the numerical suction velocity, their sensitivity is the same, meaning that in flow regimes without piping conditions, changes in the differential pressure are proportional to variations in $\alpha$.

In the experiments, the pressure difference between the inside and outside of the bucket causes a suction force on the bucket in the vertical direction. In the simulation, this force is applied to the bucket wall segment. Note that in both the simulations and the experiments, no gravitational force is acting on the bucket wall. In the experiments, the gravitational force is almost exactly balanced by the counterweight. The total suction force $F_{\text{bucket}}^{\text{suction}}$ acting on the bucket due to the pressure

<figure><img src="image_15.png" /><figcaption>**Table 6**

Experimental suction rates *Q* and corresponding simulation suction velocities $u_{\text{suction}}$ using different $\alpha$ for the cases C1-C3 with fixed bucket wall.</figcaption></figure>

<table><thead><tr><th rowspan="2">Case</th><th rowspan="2">Q (L/min)</th><th colspan="3">u<sub>suction</sub> (m/s)</th></tr><tr><th>α = 0.5</th><th>α = 0.63</th><th>α = 1.0</th></tr></thead><tbody><tr><td>C1</td><td>1.4</td><td>1.49e-3</td><td>1.87e-3</td><td>2.97e-3</td></tr><tr><td>C2</td><td>2.8</td><td>2.97e-3</td><td>3.74e-3</td><td>5.94e-3</td></tr><tr><td>C3</td><td>4.2</td><td>4.46e-3</td><td>5.61e-3</td><td>8.91e-3</td></tr></tbody></table>

difference $\Delta P$ between the inside and outside of the bucket is computed as:

$$F_{\text{bucket}}^{\text{suction}} = \Delta P \cdot A_{\text{bucket}} \quad (40)$$

Since the computational domain covers only a wall segment of the bucket, the force $F_{\text{segment}}^{\text{suction}}$ on the bucket wall segment is scaled according to the ratio of $d = 2.1$ cm to the bucket circumference of $2\pi R_{\text{bucket}} = 31.4$ cm:

$$F_{\text{segment}}^{\text{suction}} = F_{\text{bucket}}^{\text{suction}} \cdot (2.1 \text{ cm}/31.4 \text{ cm}). \quad (41)$$

### 4.3.2. Fixed bucket

For the fixed bucket condition, simulation and experimental results are compared for three cases with varying outflow velocities (C1-C3). Table 6 provides the suction rates of the experiments and the corresponding suction velocities for $\alpha = 0.5$, $\alpha = 0.63$, and $\alpha = 1.0$.

Fig. 12 illustrates the simulated temporal differential pressure evolution using $\alpha = 0.5$, $\alpha = 0.63$, and $\alpha = 1.0$, along with the experimental results for C1-C3. In C1, the differential pressure in the experiment remains constant at around 400 Pa over the time period of 15 s. The corresponding simulated differential pressure is constant over time as well for $\alpha = 0.5$, $\alpha = 0.63$, and $\alpha = 1.0$ at around 300 Pa, 400 Pa, and 600 Pa, respectively. The simulated differential pressure for $\alpha = 0.63$ matches the experiment perfectly in C1 since the calibration factor of $\alpha = 0.63$ is empirically determined such that the simulated differential pressure matches the experimental data for C1. When using $\alpha = 0.5$, the experimental data is significantly underestimated. In contrast, when using $\alpha = 1.0$, the experimental data is significantly overestimated. According to Darcy's law, the increase in differential pressure is expected to scale with the applied suction velocity. This behavior can be observed for C1, where an increase of $\alpha$ results in a corresponding increase in differential pressure. Similar to C1, the differential pressure in C2 remains constant in the experiment but is larger at around 600 Pa. The simulated differential pressure remains constant for $\alpha = 0.5$ and $\alpha = 0.63$ at around 600 Pa and 650 Pa, respectively. In contrast, a significant drop in differential pressure can be observed for $\alpha = 1.0$, which corresponds to piping conditions. The significant pressure drop indicates that choosing $\alpha = 1.0$ leads to an overestimation of piping conditions in the simulation. This observation for $\alpha = 1.0$ can be confirmed in C3, where the differential pressure drops to around 300 Pa already at $t = 5$ s. In contrast, the experimental pressure initially increases up to around 650 Pa, decreases just below 600 Pa at $t = 7$ s, followed by a steeper decrease in the upcoming two
<figure><img src="image_16.png" /><figcaption>Fig. 12. Differential pressure vs. simulation time using $\alpha = 0.5$, $\alpha = 0.63$, $\alpha = 1.0$, and the experimental data for fixed bucket conditions corresponding to (left) C1, (middle) C2, and (right) C3.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Soil states of C1-C3 in close proximity to the bucket wall at $t = 15$ s; colors according to the grain displacement accumulated over $15$ s.</figcaption></figure>

seconds down to around 400 Pa, and finally decreases slowly down to around 350 Pa at $t = 15$ s. Although the experiment shows piping, no piping condition is present when using $\alpha = 0.5$. However, for $\alpha = 0.63$, a good qualitative, as well as quasi-quantitative agreement between the numerically obtained results and the experiment can be observed. These results lead to the conclusion that $\alpha = 0.63$ yields the best agreement between the simulations and the experimentally obtained reference data. Therefore, $\alpha = 0.63$ will be used in the further scope of this paper.

In Fig. 13, a visualization of the simulated soil states of C1-C3 at $t = 15$ s is presented with colors according to the accumulated grain displacement. For the smallest suction velocity (C1), no grain displacement is visible, which confirms the choice of C1 as the quasi-static case for determining the calibration factor $\alpha = 0.63$. For the case with increased suction velocity (C2), a small heave occurs on the outflow side, and some particles next to the bucket wall are fluidized, i.e., suspended by the fluid, resulting in a small grain displacement. In C3, particles are dragged underneath the bucket wall tip, creating a channel, contributing to the significant drop in differential pressure observed in Fig. 12. The formation of such a channel, in combination with a significant drop in differential pressure, is characteristic of a fully developed piping condition. Many particles on the right next to the bucket wall are in a fluidized state.

### 4.3.3. Unrestrained bucket

For the unrestrained bucket, two simulations and corresponding experiments were conducted (C4-C5) that differ in terms of outflow velocity and suction force on the bucket. In contrast to C1-C3, the bucket position is not fixed anymore but is allowed to move in the vertical direction. Table 7 provides the suction rates and differential pressure values of the experiments and the corresponding suction velocities and suction forces for both cases. Fig. 14 compares the temporal evolution

of the bucket embedment and differential pressure, as well as the differential pressure vs. bucket embedment, between the simulations and experiments for C4-C5. The initial bucket embedment was around $e = 20$ mm. In the left pane of Fig. 14, the bucket embedment remains constant over time in C4 at around 21 mm both in the simulation and experiment. When increasing the suction velocity and suction force (C5), the bucket embedment increases with a nearly constant velocity over the time period of $15$ s both in the simulation and the experiment, finally reaching an embedment of around 40 mm and 34 mm, respectively. The resulting bucket embedment is expected to differ to some extent between the simulation and experiment for many reasons, particularly due to the differences in the investigated setups (bucket wall segment vs. entire bucket, bucket curvature is ignored, geometry of the granular sample, etc.) and initial conditions (porosity, stress state, etc.) among other factors like boundary effects. Concerning the influence of the adopted set of simulation parameters on the computed results, we generally observed a low sensitivity of the resulting bucket embedment with respect to the particle-wall friction coefficient and spring-dashpot parameters. A detailed parametric study is, however, out of the paper's scope intended here. In the middle pane of Fig. 14, the temporal evolution of the differential pressure is presented. While the differential pressure remains constant in C4 at around 400 Pa both in the simulation and experiment, it increases over time in C5 to around 1100 Pa in the simulation and 750 Pa in the experiment. A higher simulated differential pressure is expected in C5 to some extent since the bucket embedment is larger in the simulation as well, resulting in a longer path around the bucket wall segment, which is in accordance with Darcy's law. This is confirmed by the right panel of Fig. 14, where the differential pressure vs. bucket embedment is presented. This analysis compensates for the higher simulated differential pressure due to the larger bucket embedment, which makes the difference between the simulation and experiment smaller. Furthermore, an alternative simulation $C5^*$ is introduced, which differs from C5 only in the initial
<figure><img src="image_18.png" /><figcaption>Fig. 14. (left) Bucket embedment vs. time; (middle) differential pressure vs. time; (right) differential pressure vs. bucket embedment for C4-C5.</figcaption></figure>

<caption>**Table 7**</caption>

Experimental suction rates *Q* and differential pressure *ΔP*; suction velocities *u*<sub>suction</sub> using *α* = 0.63 and suction force acting on the bucket wall segment in the simulation *F*<sub>suction segment</sub> for the unrestrained cases.

<table><thead><tr><th rowspan="2">Case</th><th colspan="2">Experiment</th><th colspan="2">Simulation</th></tr><tr><th>Q (L/min)</th><th>ΔP (Pa)</th><th>u<sub>suction</sub> (m/s)</th><th>F<sub>suction segment</sub> (N)</th></tr></thead><tbody><tr><td>C4</td><td>1.6</td><td>425</td><td>2.14e-3</td><td>2.23e-1</td></tr><tr><td>C5</td><td>3.1</td><td>600</td><td>4.14e-3</td><td>3.15e-1</td></tr></tbody></table>

seabed conditions. In C5, the usual wished-in-place condition is used for the bucket, i.e., it is created to be already embedded into the soil. In contrast, in C5*, the bucket is gently pushed into the soil before the simulation starts, creating small uplifts on both sides of the bucket wall, thereby leading effectively to a slightly higher initial embedment *e* than in C5. In C5*, the differential pressure is a bit higher than in C5 for the following reasons. First, due to the higher bucket embedment. Second, the soil is initially compressed, leading to more clogging and less porosity and, therefore, higher differential pressure. Due to the low bucket embedment, the differences between C5 and C5* are relatively small. The influence of the wished in-place condition and, more broadly, of the installation effects is generally non-negligible for many geotechnical foundation systems and remains an open topic in the literature.

Fig. 15 illustrates the simulated soil states of C4-C5 at *t* = 15 s with colors according to the accumulated grain displacement. No grain displacement is visible in C4 since the suction velocity is too small, similar to C1. In contrast to C4, the bucket wall segment has significantly penetrated the soil in C5. A heave/soil plug occurs on the outflow side (Ragni et al., 2020), resulting in small grain displacements close to the bucket wall on the outflow side but no particle fluidization. However, even though the suction velocity of C5 is larger than in C2, significantly less accumulated displacement is present in C5 than in C2 (see Fig. 13), leading to the conclusion that the penetration of the bucket wall inhibits the occurrence of fluidization.

## 4.4. Micromechanical insights

This extended micromechanical analysis of the numerical results includes insights into the contact forces, particle velocities, and hydrodynamic forces for C1-C5. In Fig. 16, the simulated contact forces vs. time are plotted for C1-C5. Constant contact forces can be observed for C1, C2, and C4 since there are no significant changes in the soil states in those three cases, except that C4 has larger contact forces than C1 and C2 since the bucket wall is unrestrained and pushes in the downward direction, increasing the stress state of C4. In contrast, C3 and C5 experience a drop in contact forces, mainly caused by fluidization (C3) and soil plug expansion (C5) of the particles on the

<figure><img src="image_19.png" /><figcaption>Fig. 15. Soil states of C4-C5 in close proximity to the bucket wall at *t* = 15 s; colors according to the grain displacement accumulated over 15 s with identical quantitative limits as in Fig. 13.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Contact forces, i.e., for each particle, sum over magnitude of all contact forces with colliding particles averaged over all particles vs. time for C1-C5.</figcaption></figure>

outflow side, significantly reducing their contact forces. As for C4, the average contact force in the unrestrained case (C5) is larger than in the fixed case (C3) since the bucket wall pushes onto the particles.

Fig. 17 illustrates the contact forces at *t* = 15 s for C1-C5. For the fixed bucket and low suction velocity (C1), the contact forces mainly counteract gravity. On the outflow side, a decrease in contact forces was observed. For an increased velocity (C2), those contact forces on the outflow side are further reduced, but larger contact forces on the inflow side are visible compared to C1, caused by the increasing drag of
<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Soil states of C1-C5 at *t* = 15 s in close proximity to the bucket wall; colors indicating contact forces (magnitude); suction velocity is given as white text.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** Particle velocity vs. time for C1-C5.</figcaption></figure>

the suction. For the high-velocity fixed bucket simulation (C3), due to the piping condition and large fluidization, the contact forces become even smaller than in C1. If the bucket is unrestrained, large contact forces are present on the inflow side and especially below the bucket wall tip, caused by two effects. First, the applied suction velocity drags the particles underneath the wall tip. Second, the suction force on the bucket compresses the particles underneath the wall tip. Since both effects are more dominant in C5 than in C4, the resulting contact forces underneath the wall tip are higher as well. The decrease of the effective stresses on the outflow side and, hence, a decrease of the side friction inside the bucket is caused by the upward flow in the soil plug (Sturm, 2017).

Fig. 18 presents the average particle velocity vs. time for C1-C5. Either no or very small particle velocities are present for the cases with small particle displacement, i.e., C1, C2, and C4. In contrast, there are higher particle velocities in C3 and C5. An initial spike can be observed in C3 due to the incipient of piping, followed by a decrease and convergence due to the fluidization on the outflow side. There is little piping and fluidization present in C5. Therefore, the average particle velocity is smaller than in C3 and mainly originates from the soil plug expansion. This aligns with the accumulated grain displacement, which is significantly lower in C5 than in C3 (see Fig. 13 and Fig. 15).

Fig. 19 illustrates the time-averaged velocity vs. the coordination number for individual particles of C3. An inverse correlation between the coordination number and the particle velocity is observed, i.e., particles with higher velocities tend to have a lower coordination number and vice versa. This inverse correlation is expected since particles with higher velocities are, on average, in contact with fewer neighboring particles, especially once they are fluidized, and thus a smaller coordination number.

Fig. 20 shows the average hydrodynamic force on the particles vs. time and vs. differential pressure for C1-C5. The development of the

<figure><img src="image_23.png" /><figcaption>**Fig. 19.** Scatter plot with time-averaged particle velocity vs. coordination number *z* for individual particles of C3. Each point of the scatter plot represents one particle, but only every 10th particle is visualized.</figcaption></figure>

hydrodynamic forces highly correlates with the differential pressure presented in Fig. 12 and Fig. 14. The overestimation is confirmed when plotting the hydrodynamic forces vs. the differential pressure, which leads to a nearly perfect linear relationship for C1-C4, in agreement with the global force balance in the porous medium. In contrast, in C5, there is a deviation from the linear relationship between hydrodynamic force and the differential pressure between 800 Pa and 1100 Pa, which is precisely the range in which the bucket wall penetrates the soil. The non-linear relation of hydrodynamic forces and differential pressure is caused by the increase of bucket embedment in C5. This increase changes the length of the flow around the bucket wall segment through the porous media, but the linear relationship only holds for a fixed length.

Finally, Fig. 21 visualizes the hydrodynamic forces acting on the particles surrounding the bucket wall at the onset of piping in C3. The largest hydrodynamic forces can be observed close to the bucket wall where the fluid velocity is the highest, especially on the inflow side and underneath the bucket wall tip, dragging the particles in the flow direction. Nearly no hydrodynamic forces are present on the outflow side due to the fluidization of the particles beginning.

## 5. Conclusion

In this paper, a comprehensive three-dimensional, fully-resolved coupled model for a local micro-mechanical analysis of the installation of a suction bucket foundation was introduced. The LBM was employed for the Eulerian fluid flow, whereas the Lagrangian particles were modeled using the DEM. The methodology was validated using three well-established benchmarks, namely a fixed sphere surrounded by Stokes flow, a fixed sphere exposed to shear flow, and a freely settling
<figure><img src="image_24.png" /><figcaption>Fig. 20. (left) Hydrodynamic force vs. time; (right) hydrodynamic force vs. differential pressure for C1-C5.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 21. (left) Hydrodynamic forces on particles in close proximity to the bucket wall at the critical time step $t = 1.15$ s, i.e., the onset of the pressure drop, in C3; (right) same state visualized by arrows.</figcaption></figure>

sphere under gravity. A subsequent cross-comparison of the presented numerical model with an experimental study of suction-assisted bucket settling in granular soil was presented, both for a fixed and unrestrained bucket. A calibration factor $\alpha$ was introduced, which accounts for the geometrical differences between the simulation domain and the experimental setup. The differential pressure, a key indicator for piping erosion, and the bucket embedment vs. time were compared between the numerical and experimental results. Finally, an in-depth analysis of the micromechanics was provided, including insights into the contact forces, particle velocities, and hydrodynamic forces.

The presented results show that:

* It is computationally feasible to conduct three-dimensional fully-resolved micro-mechanical simulations of piping erosion in relevant scenarios at real scale using the LBM and DEM.

* With little macromechanical model assumptions, the computed results are not only physically sound but also match experimental data, in good agreement with theoretical expectations from the suction bucket literature.

* The onset of granular fluidization (piping erosion) experimentally observed in simple setups with a fixed bucket can be well reproduced numerically, showing a good quantitative match of

the temporal evolution of suction pressure with the experimental data.

* For the more complex scenarios with an unrestrained bucket, the computed evolution of bucket embedment and differential pressure still shows a very good qualitative agreement with the experimental data and remains within reasonable quantitative margins given the obvious, generally unavoidable, modeling simplifications.

* The mechanical resistance during the suction penetration is concentrated outside the bucket and underneath the bucket wall tip, whereas the plug inside the bucket is loosened, leading to a decrease in contact forces and friction.

Within the limited scope of the presented study, it is possible to derive the following lessons:

* Given the unavoidable and significant differences between the assumed and the real conditions, a careful consideration of experimental data seems essential for calibrating the model.

* For a general quantitative match with the available experimental data, it appears pertinent to first estimate an initial calibration factor out of a quasi-static reference case, such that the computed differential pressure matches the experimental one.
* The results show a higher sensitivity to the physical installation conditions, e.g., suction velocity and initial bucket embedment, than to the particular choice of the DEM parameters.

* In practical terms of computational effort, it appears desirable to use at least one supercomputing node, e.g., eight GPUs for about 20 h, for the simulation of setups comparable to those presented here.

Furthermore, our future work will also envisage significantly larger simulation setups using large-scale supercomputers with hundreds of nodes, allowing the consideration of suction bucket prototypes in realistic scales in the context of offshore wind turbines. The presented modeling approach provides access and detailed insights into the physics of fluid-saturated granular systems and, in this case, the triggering conditions of piping erosion. The proposed model demonstrates very promising capabilities for analyzing geotechnical problems such as piping erosion. In general, the full-scale simulation of entire geotechnical structures is, in principle, out of reach due to the enormous computational requirements with the fully-resolved discrete approach. Nevertheless, local phenomena with moderate spatial extents taking place in representative or dominant parts of the foundation may well be analyzed in real size. In this sense, other similar potential applications involving granular fluidization could include, for instance, local scenarios of backward erosion during the onset of a dam or levee breach, soil liquefaction, e.g., around dynamically driven displacement piles, or hydraulic failure problems in excavation pits below the phreatic level in cohesionless soils.

### CRediT authorship contribution statement

**Samuel Kemmler:** Writing – original draft, Visualization, Validation, Software, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Pablo Cuéllar:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. **Antoni Artinov:** Writing – review & editing, Conceptualization. **Li-Hua Luu:** Writing – review & editing, Conceptualization. **Abbas Farhat:** Writing – review & editing, Data curation. **Pierre Philippe:** Writing – review & editing, Conceptualization. **Christoph Rettinger:** Writing – review & editing, Supervision, Resources, Funding acquisition.

### 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.

### Acknowledgments

This work has received funding from the European High Performance Computing Joint Undertaking (JU) and Sweden, Germany, Spain, Greece, and Denmark under grant agreement No 101093393. The funding source had no direct involvement in this paper.

The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU). The hardware is funded by the German Research Foundation (DFG).

The authors gratefully acknowledge Susanne Malheiros, Uwe Wössner, and the HLRS for creating the visualization of the graphical abstract from the simulation data provided to them using the scientific visualization software Vistle (Aumüller, 2019).

The authors gratefully acknowledge the Norwegian Geotechnical Institute (NGI) and Ørsted for providing images of the Borkum Riffgrund 2 installation.

### Data availability

Code and data are available on Zenodo: https://zenodo.org/records/14170253.

### References

Aumüller, M., 2019. Hybrid remote visualization in immersive virtual environments with vistle. In: Eurographics Symposium on Parallel Graphics and Visualization. The Eurographics Association, http://dx.doi.org/10.2312/PGV.20191113.

Bauer, M., Eibl, S., Godenschwager, C., Kohl, N., Kuron, M., Rettinger, C., Schornbaum, F., Schwarzmeier, C., Thönnes, D., Köstler, H., Rüde, U., 2021. Walberla: A block-structured high-performance framework for multiphysics simulations. Comput. Math. Appl. 81, 478–501. http://dx.doi.org/10.1016/j.camwa.2020.01.007, arXiv:1909.13772.

Benseghier, Z., Cuéllar, P., Luu, L.H., Bonelli, S., Philippe, P., 2020. A parallel GPU-based computational framework for the micromechanical analysis of geotechnical and erosion problems. Comput. Geotech. 120, 103404. http://dx.doi.org/10.1016/j.compgeo.2019.103404.

Biegert, E., Vowinckel, B., Meiburg, E., 2017. A collision model for grain-resolving simulations of flows over dense, mobile, polydisperse granular sediment beds. J. Comput. Phys. 340, 105–127. http://dx.doi.org/10.1016/j.jcp.2017.03.035.

Cerfontaine, B., Ciantia, M., Brown, M., White, D., Sharif, Y., 2023. DEM study of particle scale effect on plain and rotary jacked pile behaviour in granular materials. Comput. Geotech. 161, 105559. http://dx.doi.org/10.1016/j.compgeo.2023.105559.

Ciantia, M.O., Arroyo, M., Butlanska, J., Gens, A., 2016. DEM modelling of cone penetration tests in a double-porosity crushable granular material. Comput. Geotech. 73, 109–127. http://dx.doi.org/10.1016/j.compgeo.2015.12.001.

Cotter, O., 2009. Installation of Suction Caisson Foundations for Offshore Renewable Energy Structures (Ph.D. thesis). University of Oxford.

Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Géotechnique 29 (1), 47–65. http://dx.doi.org/10.1680/geot.1979.29.1.47.

Delenne, J.Y., EI Youssoufi, M.S., Cherblanc, F., Bénet, J.C., 2004. Mechanical behaviour and failure of cohesive granular materials. Int. J. Numer. Anal. Methods Geomech. 28 (15), 1577–1594. http://dx.doi.org/10.1002/nag.401.

Froio, F., Callari, C., Rotunno, A.F., 2019. A numerical experiment of backward erosion piping: Kinematics and micromechanics. Meccanica 54 (14), 2099–2117. http://dx.doi.org/10.1007/s11012-019-01071-7.

Fukumoto, Y., Ohtsuka, S., 2018. 3-d direct numerical model for failure of non-cohesive granular soils with upward seepage flow. Comput. Part. Mech. 5 (4), 443–454. http://dx.doi.org/10.1007/s40571-017-0180-5.

Fukumoto, Y., Yang, H., Hosoyamada, T., Ohtsuka, S., 2021. 2-d coupled fluid-particle numerical analysis of seepage failure of saturated granular soils around an embedded sheet pile with no macroscopic assumptions. Comput. Geotech. 136, 104234. http://dx.doi.org/10.1016/j.compgeo.2021.104234.

Ginzburg, I., Verhaeghe, F., 2008. Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions. Commun. Comput. Phys.

Holzer, M., Bauer, M., Köstler, H., Rüde, U., 2021. Highly efficient lattice Boltzmann multiphase simulations of immiscible fluids at high-density ratios on CPUs and GPUs through code generation. Int. J. High Perform. Comput. Appl. 35 (4), 413–427. http://dx.doi.org/10.1177/10943420211016525.

Jones, B.D., Williams, J.R., 2017. Fast computation of accurate sphere-cube intersection volume. Eng. Comput. 34 (4), 1204–1216. http://dx.doi.org/10.1108/EC-02-2016-0052.

Kay, S., Gourvenec, S., Palix, E., Alderlieste, E., 2021. Intermediate Offshore Foundations, first ed. CRC Press, London, http://dx.doi.org/10.1201/9780429423840.

Kemmler, S., Cuéllar, P., Rettinger, C., Köstler, H., 2025a. A fluid-solid coupled micromechanical simulation for the analysis of piping erosion during the seabed installation of a suction bucket foundation. IOP Conf. Ser.: Earth Environ. Sci. 1480 (1), 012024. http://dx.doi.org/10.1088/1755-1315/1480/1/012024.

Kemmler, S., Rettinger, C., Rüde, U., Cuéllar, P., Köstler, H., 2025b. Efficiency and scalability of fully-resolved fluid-particle simulations on heterogeneous CPU-GPU architectures. Int. J. High Perform. Comput. Appl. 39 (3), 345–363. http://dx.doi.org/10.1177/10943420241313385.

Khirevich, S., Ginzburg, I., Tallarek, U., 2015. Coarse- and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings. J. Comput. Phys. 281, 708–742. http://dx.doi.org/10.1016/j.jcp.2014.10.038.

Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., Viggen, E.M., 2017. The lattice Boltzmann method: principles and practice. In: Graduate Texts in Physics, Springer International Publishing, Cham, http://dx.doi.org/10.1007/978-3-319-44649-3.

Ladd, A.J.C., Verberg, R., 2001. Lattice-Boltzmann simulations of particle-fluid suspensions. J. Stat. Phys. 104 (5), 1191–1251. http://dx.doi.org/10.1023/A:1010414013942.
Luding, S., 2008. Cohesive, frictional powders: Contact models for tension. Granul. Matter 10 (4), 235-246. http://dx.doi.org/10.1007/s10035-008-0099-x.

Noble, D.R., Torczynski, J.R., 1998. A lattice-Boltzmann method for partially saturated computational cells. Internat. J. Modern Phys. C 09 (08), 1189-1201. http://dx.doi.org/10.1142/S0129183198001084.

Penner, N., Grießmann, T., Rolfes, R., 2020. Monitoring of suction bucket jackets for offshore wind turbines: Dynamic load bearing behaviour and modelling. Mar. Struct. 72, 102745. http://dx.doi.org/10.1016/j.marstruc.2020.102745.

Potyondy, D., Cundall, P., 2004. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 41 (8), 1329-1364. http://dx.doi.org/10.1016/j.ijrmms.2004.09.011.

Ragni, R., Bienen, B., Stanier, S., O'Loughlin, C., Cassidy, M., 2020. Observations during suction bucket installation in sand. Int. J. Phys. Model. Geotech. 20 (3), 132-149. http://dx.doi.org/10.1680/jphmg.18.00071.

Rettinger, C., 2023. Fully Resolved Simulation of Particulate Flows with a Par- allel Coupled Lattice Boltzmann and Discrete Element Method (Ph.D. thesis). Friedrich-Alexander-Universität Erlangen-Nürnberg.

Rettinger, C., Rüde, U., 2017. A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann simulations. Comput. & Fluids 154, 74-89. http://dx.doi.org/10.1016/j.compfluid.2017.05.033.

Rettinger, C., Rüde, U., 2018. A coupled lattice Boltzmann method and discrete element method for discrete particle simulations of particulate flows. Comput. & Fluids 172, 706-719. http://dx.doi.org/10.1016/j.compfluid.2018.01.023.

Rettinger, C., Rüde, U., 2022. An efficient four-way coupled lattice Boltzmann - discrete element method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 453, 110942. http://dx.doi.org/10.1016/j.jcp.2022.110942.

Sangani, A., Acrivos, A., 1982. Slow flow through a periodic array of spheres. Int. J. Multiph. Flow 8 (4), 343-360. http://dx.doi.org/10.1016/0301-9322(82)90047-7.

Schwarzmeier, C., Rettinger, C., Kemmler, S., Plewinski, J., Núñez-González, F., Köstler, H., Rüde, U., Vowinckel, B., 2023. Particle-resolved simulation of antidunes in free-surface flows. J. Fluid Mech. 961, R1. http://dx.doi.org/10.1017/jfm.2023.262.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations. Mon. Weather Rev. 91 (3), http://dx.doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2.

Sturm, H., 2017. Design aspects of suction caissons for offshore wind turbine foundations. In: Unearth the Future, Connect beyond. pp. 45-63.

Swope, W.C., Andersen, H.C., Berens, P.H., Wilson, K.R., 1982. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: application to small water clusters. J. Chem. Phys. 76 (1), 637-649. http://dx.doi.org/10.1063/1.442716.

Tang, H., Song, R., Dong, Y., Song, X., 2019. Measurement of restitution and friction coefficients for granular particles and discrete element simulation for the tests of glass beads. Materials 12 (19), 3170. http://dx.doi.org/10.3390/ma12193170.

Ten Cate, A., Nieuwstad, C.H., Derksen, J.J., Van Den Akker, H.E.A., 2002. Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity. Phys. Fluids 14 (11), 4012-4025. http://dx.doi.org/10.1063/1.1512918.

Tran, D., Prime, N., Froio, F., Callari, C., Vincens, E., 2017. Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling. Eur. J. Environ. Civ. Eng. 21 (7-8), 960-987. http://dx.doi.org/10.1080/19648189.2016.1248794.

Wang, D., Leonardi, C.R., Aminossadati, S.M., 2018. Improved coupling of time integration and hydrodynamic interaction in particle suspensions using the lattice Boltzmann and discrete element methods. Comput. Math. Appl. 75 (7), 2593-2606. http://dx.doi.org/10.1016/j.camwa.2018.01.002.

Yang, Y., Bao, J., Yang, H., Pang, Z., Wang, Y., 2023. Microscopic analysis of the influence of soil properties on the suction bucket installation in sand based on the CFD-DEM model. Comput. Geotech. 156, 105249. http://dx.doi.org/10.1016/j.compgeo.2023.105249.

Yu, H., Girimaji, S.S., Luo, L.S., 2005. DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method. J. Comput. Phys. 209 (2), 599-616. http://dx.doi.org/10.1016/j.jcp.2005.03.022.

Zeng, L., Najjar, F., Balachandar, S., Fischer, P., 2009. Forces on a finite-sized particle located close to a wall in a linear shear flow. Phys. Fluids 21 (3), 033302. http://dx.doi.org/10.1063/1.3082232.

Zhang, P., Mu, L., Lu, Y., Huang, M., Sun, J., 2024. Microscopic insights into suction bucket installation in sand: coupled coarse-grained CFD-dem simulations. Comput. Geotech. 167, 106060. http://dx.doi.org/10.1016/j.compgeo.2023.106060.