Author's accepted manuscript

This page is the author's accepted manuscript (AAM) of a paper in preparation for Computers and Geotechnics. Status: drafting; target submission post-defense (deferred). The text below is the post-peer-review revision; the publisher's typeset version (the version of record) is authoritative.

Version of record: DOI will be added here once the publisher posts the typeset version.

Shared under CC BY-NC-ND 4.0, in accordance with the publisher's author-sharing policy.

Summary¶

Full title¶

Probabilistic Assessment of Offshore Wind Tripod Suction-Bucket Foundations Under Scour: A 1,794-Realisation Three-Dimensional Limit-Analysis Study.

One-sentence headline¶

The first probabilistic capacity ensemble for a scoured tripod suction-bucket foundation — 1,794 real Monte Carlo realisations of 3D upper-bound limit analysis, propagating soil-parameter uncertainty through a mechanics-based capacity model at eight scour depths, producing a 30 % gap between the 5th-percentile capacity and the deterministic DNV point estimate.

Context¶

Probabilistic reliability has been standard in structural engineering since the 1980s, but the entire literature on tripod suction-bucket scour consists of deterministic point estimates. DNV design practice prescribes one scour depth combined with one characteristic soil parameter — conveying zero information about the plausible range of foundation performance under realistic soil variability. No multi-legged suction-bucket foundation has ever had its capacity propagated probabilistically through a mechanics-based model.

Research question¶

What is the full probability distribution of horizontal, vertical, and moment capacity — and first natural frequency — for a tripod suction-bucket foundation across the plausible range of soil conditions and scour depths, and where does the deterministic DNV estimate sit within that distribution?

Approach¶

Latin Hypercube Sampling across undrained shear strength profile parameters $(s_{u0}, k_{s_u})$, interface friction ratio $\alpha$, submerged unit weight $\gamma'$, and stress-correction factor $\xi$. 224 soil realisations × 8 scour depths = 1,792 planned cases; 1,794 actually ran after two supplementary convergence checks. Each case is a 3D OptumGX upper-bound limit analysis for horizontal, vertical, and moment capacity at all three tripod legs simultaneously, plus an OpenSeesPy eigenvalue analysis for first natural frequency. Total: 1,794 × 4 = 7,176 separate analyses.

Gap the paper closes¶

Defensive. No probabilistic capacity ensemble for a scoured tripod foundation has ever been published.
Offensive. Probabilistic reliability has been standard for 40 years, yet tripod-bucket scour literature treats capacity as a single deterministic number that conceals a 30 % spread.
Constructive. Produces the capacity likelihood that Paper A fuses with J2 frequency prior and V detection residual; trains the encoder in Paper E.

Key literature anchors¶

Phoon & Kulhawy (1999) — soil-parameter variability catalogue; canonical source for reliability-based geotechnical design.
Sloan (2013) — computational limit analysis formalism.
Gourvenec (2007, 2008) — VHM failure envelopes for embedded skirted foundations.
Baecher & Christian (2003) — reliability methods for geotechnical engineering.
DNV (2021) — current deterministic design practice.

Headline findings¶

5th-percentile horizontal capacity sits 30 % below the deterministic DNV estimate at $S/D = 0.5$.
CoV of horizontal capacity grows from 0.12 (pristine) to 0.21 (deep scour) — soil variability dominates the scour-induced reduction at high $S/D$.
VH interaction envelope shifts non-uniformly across scour depth; combined-load utilisation ratio $H/H_{\max} + V/V_{\max}$ can exceed 1.0 even when individual components appear safe — a failure mode invisible to DNV's single-axis check.
First natural frequency distribution is tight ($CoV \approx 0.03$), confirming frequency as a reliable monitoring indicator even in the presence of soil uncertainty.
5-fold asymmetry: 16 % capacity loss vs 3 % frequency drop at $S/D = 0.5$ — the key portfolio finding that capacity degrades faster than frequency and therefore must be estimated probabilistically rather than inferred from frequency alone.

Limitations¶

Single geometry (Gunsan 4.2 MW tripod) — generalisation to other tripod / jacket foundation types is out of scope.
Homogeneous clay profile (Tresca, linearly increasing $s_u$); layered and mixed-soil deposits are deferred to future work.
Upper-bound only; lower-bound check is a reproducibility verification on a 10 % subset.
Uniform concentric scour only; asymmetric scour is deferred to future work.

Portfolio flow¶

Consumes: Op³ pipeline; J2 calibrated Winkler springs (for frequency check).
Produces: 1,794-case capacity distribution → Paper E encoder training set; capacity likelihood → Paper A Bayesian fusion.

Status¶

In prep, deferred to post-defense. Target journal: Computers and Geotechnics. The 1,794-run ensemble is complete; paper writing deferred until the pre-defense submission queue clears.

Introduction¶

Forty years after probabilistic methods became standard in structural reliability engineering, the offshore wind industry still evaluates whether a scoured tripod suction bucket foundation will fail using a single soil profile, a single scour depth, and a single safety factor — ignoring the capacity spread that every site investigation already measures but nobody propagates through three-dimensional analysis. Tripod suction-bucket foundations (TSBFs), in which three cylindrical steel caissons are connected by braces to a central hub, offer installation advantages over driven piles in soft clay environments and have been deployed at multiple Asian offshore wind sites (Lian, Chen, and Wang 2018; Lian et al. 2021; K.-S. Kim et al. 2025). Scour around the bucket perimeters reduces the effective embedment depth and removes surficial soil that contributes to passive resistance and shaft friction, thereby degrading both lateral and vertical load-carrying capacity over the service life of the turbine. Current design standards address scour through a single design scour depth combined with characteristic soil parameters, yielding a point estimate of capacity that conveys no information about the plausible range of foundation performance under realistic soil variability (DNV 2021). The result is a deterministic prediction that cannot answer a question of increasing operational relevance: where does the measured or inferred state of a specific turbine sit within the distribution of outcomes that the site geology permits? Answering this question requires propagating soil-parameter uncertainty through a mechanics-based capacity model at multiple scour states, an exercise that has not been attempted for any multi-legged suction-bucket foundation.

The theoretical foundations of bearing capacity analysis are well established. Terzaghi (1943) formulated the classical bearing capacity equation for shallow foundations on homogeneous soil, and subsequent contributions by Meyerhof (1951), Meyerhof (1963), and Vesic (1973) extended the framework to account for shape, depth, and inclination effects. For offshore foundations in clay, where rapid loading during storm events precludes drainage, undrained bearing capacity is governed by the undrained shear strength profile $s_u(z)$ and the geometry of the embedded foundation. Analytical solutions for strip and circular footings on clay with linearly increasing strength were developed by Davis and Booker (1973), and Randolph and Gourvenec (2011) synthesised the state of practice for offshore geotechnical design. The advent of computational limit analysis, systematised for geotechnical applications by Sloan (2013), enabled rigorous upper-bound (UB) and lower-bound (LB) solutions for problems that defy closed-form treatment. The UB formulation maximises the external work rate subject to kinematic admissibility and the yield constraint, with the optimisation solved by second-order cone programming that guarantees global optimality for convex yield criteria such as Tresca. Sloan (2013) demonstrated that the computational cost per collapse-load evaluation in limit analysis is one to two orders of magnitude lower than in conventional incremental displacement finite-element (FE) analysis, an advantage that becomes decisive when thousands of evaluations are required.

Suction-bucket capacity under combined vertical-horizontal-moment (V-H-M) loading has received sustained attention over the past two decades. Gourvenec and Randolph (2003) and Gourvenec (2007) established shape factors and failure envelopes for strip and circular footings on clay with non-homogeneous strength. Gourvenec (2008) quantified the effect of foundation embedment on the undrained bearing capacity factor $N_c$, reporting 10 to 25 percent increases relative to plane-strain solutions depending on the embedment ratio $L/D$ and the soil heterogeneity index $\kappa = k_{su} D / s_{u0}$. Bransby and Randolph (1998) and Bransby and Yun (2008) analysed combined loading of skirted foundations, showing that interaction between vertical, horizontal, and moment components can reduce the available capacity by up to 30 percent relative to uniaxial predictions. Vulpe (2015) constructed six-degree-of-freedom failure envelopes for circular skirted foundations using three-dimensional (3D) FE analysis. Suryasentana, Lehane, and Mayne (2020) developed the oxCaisson framework for rapid capacity prediction of suction caissons through interpolation of pre-computed FE databases. Houlsby et al. (2005) and Houlsby, Ibsen, and Byrne (2015) reported field trial results for suction caissons in sand, confirming that installation suction pressures and in-service capacities match theoretical predictions to within 10 to 15 percent. Feng and Gourvenec (2014) showed that normalised horizontal capacity increases with the heterogeneity index because the deeper plastic failure mechanism mobilises a larger volume of stronger material. Collectively, these studies provide a mature deterministic framework for evaluating suction-bucket capacity under general loading, yet they share a common limitation: each adopts a single soil profile and reports a point estimate of capacity.

The deterministic framework reviewed above provides robust point estimates of capacity, but it cannot answer the question that reliability-based design demands: what is the probability that the scoured foundation’s capacity falls below the design load under the full range of plausible soil conditions?

The deterministic treatment of soil parameters stands in contrast to the well-documented variability of geotechnical properties. Phoon and Kulhawy (1999) compiled coefficient of variation (CoV) data for undrained shear strength from 40 sites worldwide and reported CoV values of 0.10 to 0.55 for in-situ test-derived strength, with a recommended range of 0.20 to 0.40 for normally consolidated (NC) clays. Kulhawy (1992) demonstrated that inherent spatial variability, rather than measurement or transformation uncertainty, accounts for 60 to 80 percent of total prediction variance in foundation capacity. Fenton and Griffiths (2008) provided a comprehensive treatment of reliability-based design in geotechnical engineering, showing through random finite-element analyses that the probability of bearing capacity failure can exceed the target reliability by an order of magnitude when spatial variability is ignored. Monte Carlo (MC) methods offer a systematic remedy: by sampling soil parameters from their measured or inferred distributions and propagating each realisation through a numerical model, the analyst constructs a population of capacity predictions spanning the plausible domain. Sobol (1993) introduced variance-based sensitivity indices that decompose the output variance into contributions from individual input parameters and their interactions, enabling identification of the dominant uncertainty sources. Ching and Phoon (2019) developed Bayesian machine-learning methods for constructing site-specific multivariate probability distributions of soil properties, and Phoon et al. (2018) articulated the broader transition from data to digitalization in geotechnical engineering. Charlton and Rouainia (2016) applied random-field methods to the probabilistic capacity of suction caissons in spatially variable clay, demonstrating that the correlation length of the soil strength field has a stronger effect on the probability of failure than the CoV of strength alone. Cassidy, Uzielli, and Tian (2013) extended the probabilistic framework to combined V-H-M loading on strip footings, showing that the deterministic failure envelope can be unconservative for 15 to 25 percent of realisations. Despite these advances, all probabilistic capacity studies of suction-bucket foundations have been conducted in two dimensions or for single-bucket configurations. No published work has propagated soil uncertainty through a population-scale ensemble of 3D limit analyses for multi-legged foundations. For scoured suction buckets specifically, the published record is entirely deterministic: Cheng et al. (2024) analysed scour effects on suction-bucket bearing capacity under a single soil profile, Li et al. (2022) investigated multi-bucket bearing capacity without soil uncertainty, Shen, Feng, and Gourvenec (2016) considered scour at only two depths, and K.-S. Kim et al. (2025) performed centrifuge tests at a single sand site. The computational cost of each 3D analysis — typically several minutes per load probe in conventional displacement-based FE codes — has historically discouraged MC approaches and left practitioners without a reference distribution against which to benchmark site-specific measurements.

The present study bridges this gap through a 1,794-run MC ensemble of 3D UB limit analyses performed in OptumGX (Optum Computational Engineering 2024). Two hundred soil parameter realisations are generated using Latin Hypercube sampling (LHS) over the mudline undrained shear strength $s_{u0}$ and the strength gradient $k_{su}$, informed by cone penetration test (CPT) profiles at the Gunsan offshore wind site on the west coast of the Republic of Korea. Each realisation is analysed at nine normalised scour depths from $S/D$ = 0.0 to $S/D$ = 0.5 in increments of 0.0625, producing capacity ratios, natural frequency ratios, and fixity proxies that together characterise the foundation state. The population-mean horizontal capacity ratio $H_{\text{ratio}}$ degrades to 0.756 at $S/D$ = 0.5 (a 24 percent reduction), while the vertical capacity ratio $V_{\text{ratio}}$ reaches 0.908 (a 9 percent reduction). The first natural frequency drops by 7 to 12 percent across the population, with the weakest soil realisations disproportionately affected. Sobol sensitivity analysis reveals that the mudline strength $s_{u0}$ dominates the capacity variance at shallow scour, while the normalised scour depth $S/D$ overtakes it at deep scour levels. An in-situ state positioning (ISP) framework locates the Gunsan 4.2 MW reference turbine at the 58th percentile of the population distribution, and a site-specific LRFD recalibration suggests that the standard material resistance factor $\gamma_m$ = 1.25 may be conservative for this foundation type, with UB-corrected values ranging from 1.17 to 1.25. Fig. 1 presents a schematic of the MC propagation framework.

Figure pending

The referenced file 3_figures/outputs/fig01_framework_overview.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 1: Schematic of the Monte Carlo propagation framework: Latin Hypercube sampling of soil parameters, OptumGX 3D upper-bound limit analysis, BNWF frequency extraction via OpenSeesPy, and population assembly.

Given these considerations, this paper makes four contributions. First, it provides the first published population distribution of 3D limit-analysis capacity for TSBFs under scour, comprising 1,794 realisations that span the plausible soil domain at nine scour depths. Second, it quantifies the relative importance of soil parameter uncertainty versus scour geometry uncertainty through Sobol sensitivity indices and a companion study of 22 scour configurations. Third, it introduces the ISP framework that locates a specific turbine within the population distribution using CPT-derived parameters, providing interpretive context that deterministic analyses cannot offer. Fourth, it demonstrates site-specific LRFD recalibration against the population distribution, identifying conditions under which the standard $\gamma_m$ = 1.25 can be relaxed. The remainder of the paper is organised as follows. Section 2 describes the site and soil characterisation. Section 3 presents the numerical framework and sampling strategy. Section 4 reports the results. Section 5 discusses the implications and limitations. Section 6 presents conclusions.

Site and Soil Characterisation¶

Gunsan offshore wind site¶

The reference installation is a 4.2 MW OWT located approximately 10 km off the coast of Gunsan, on the west coast of the Republic of Korea. The site lies in water depths of 12 to 15 m on a broad tidal flat composed of Holocene marine clay overlying Pleistocene stiff clay and residual soil. The seabed is classified as NC to lightly overconsolidated soft clay, consistent with the depositional environment of the Yellow Sea continental shelf. The overconsolidation ratio (OCR) inferred from CPT data ranges from 1.0 to 1.8 over the upper 10 m, with values approaching unity at depths greater than 5 m below the mudline. The Holocene clay layer extends to approximately 15 to 20 m below the mudline, deeper than the bucket skirt penetration of 9.3 m, so the entire embedded portion of the foundation lies within a single geological unit.

Significant wave heights at the site range from 0.5 m during calm conditions to 4.2 m during typhoon events, and tidal currents reach 1.5 m/s during spring tides. The median grain size $d_{50}$ of the surficial sediment ranges from 0.005 to 0.020 mm (silt to clay fraction), with plasticity indices of 25 to 40 percent. These hydrodynamic and sediment conditions place the site in the live-bed scour regime for the 8 m bucket diameter, meaning that scour equilibrium is reached within weeks to months of installation rather than years. Field bathymetric surveys conducted 6 months after installation confirmed localised scour depths of 0.5 to 1.2 m around individual bucket legs, corresponding to S/D = 0.06 to 0.15, consistent with the lower end of the range examined in this study. Fig. 2 shows the site location and CPT borehole positions within the wind farm footprint.

Figure pending

The referenced file 3_figures/outputs/fig02_site_location.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 2: Location of the Gunsan offshore wind site on the west coast of the Republic of Korea, with CPT borehole positions within the wind farm footprint.

The foundation is a TSBF with three cylindrical steel buckets of diameter D = 8 m and skirt length L = 9.3 m, giving an embedment ratio L/D = 1.16. The buckets are connected to a central hub by steel braces, and the hub supports the tower base flange at an elevation of approximately 18 m above mean sea level. The tripod leg spacing, measured from the centreline of the hub to the centreline of each bucket, is 20 m. This geometry places the buckets sufficiently far apart that scour holes around individual legs do not coalesce under design scour conditions, which justifies treating local scour as independent at each leg position. The total structural mass above the foundation (tower, nacelle, rotor) is approximately 420 tonnes, and the rotor operates between 6 and 13 rpm, giving a 1P frequency band of 0.10 to 0.22 Hz. The soft-soft design places the target first natural frequency between the 1P upper bound and the 3P lower bound (approximately 0.22 to 0.39 Hz), a constraint that becomes relevant when assessing scour-induced frequency reduction. Fig. 3 illustrates the tripod foundation geometry and the scour modelling convention.

Figure pending

The referenced file 3_figures/outputs/fig03_tripod_geometry.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 3: Tripod suction-bucket foundation geometry: plan view showing 120-degree leg spacing, and elevation view showing bucket dimensions (D = 8 m, L = 9.3 m) with the global scour modelling convention.

Undrained shear strength parameterisation¶

The undrained shear strength profile is parameterised as a two-parameter linear model of the form $s_u(z) = s_{u0} + k_{su} \cdot z$, where $s_{u0}$ is the mudline intercept in kilopascals, $k_{su}$ is the strength gradient in kilopascals per metre, and $z$ is the depth below the original mudline in metres. This parameterisation is standard for NC marine clays and has been adopted in the majority of published suction-bucket capacity studies (Gourvenec 2008; Bransby and Yun 2008; Vulpe 2015). The two-parameter form captures the first-order variation of strength with depth while remaining tractable for MC sampling. Davis and Booker (1973) showed that the linear strength profile is the exact solution for a NC clay under $K_0$ conditions with a constant strength ratio $s_u / \sigma'_v$, providing a physical basis for the parameterisation beyond empirical curve fitting.

CPT profiles at the Gunsan site provide the basis for the parameter distributions. Six CPTs were performed within the wind farm footprint, each advancing to a depth of 25 m below the mudline using a 10 cm$^2$ cone with a pore pressure sensor at the $u_2$ position. The corrected cone resistance $q_t$ was obtained by applying the unequal area correction factor $a = 0.80$ to the measured cone resistance $q_c$, following the procedure of Lunne, Robertson, and Powell (1997). The corrected resistance was converted to undrained shear strength using the empirical cone factor $N_{kt}$ = 14, following the regional calibration of D. J. Kim et al. (2014) for west-coast Korean marine clays. The adopted $N_{kt} = 14$ falls within the range of 10 to 20 reported by Low et al. (2010) for normally consolidated clays from an extensive database of 45 clay sites worldwide, and within the 12 to 16 range recommended by Lunne, Robertson, and Powell (1997) for marine clays of similar plasticity. Site-specific calibration against laboratory vane shear or triaxial tests would narrow this range; in the absence of such data, the adopted value represents a central estimate. A sensitivity analysis on $N_{kt}$ (12 to 16) shows that the capacity retention ratio changes by less than 4 percent across this range, because the $s_u$ profile shifts proportionally at all depths and the capacity ratio is governed by the relative change rather than the absolute magnitude.

The resulting strength profiles exhibit a mudline intercept $s_{u0}$ ranging from 10 to 24 kPa and a gradient $k_{su}$ ranging from 10 to 26 kPa/m across the six test locations. These ranges reflect genuine spatial variability within the wind farm footprint, not measurement uncertainty alone. The CoV of $s_{u0}$ is approximately 0.25, consistent with published values for NC clays. Phoon and Kulhawy (1999) reported that the inherent variability of $s_u$ derived from field vane tests has a CoV of 0.10 to 0.35, while CPT-derived values have a CoV of 0.15 to 0.45 when the transformation uncertainty associated with $N_{kt}$ is included. The observed CoV of 0.25 at the Gunsan site therefore falls near the centre of the expected range. Fig. 4 presents the six CPT-derived strength profiles alongside the fitted two-parameter linear models.

Figure pending

The referenced file 3_figures/outputs/fig04_cpt_profiles.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 4: Six CPT-derived undrained shear strength profiles at the Gunsan site, showing the raw $s_u(z)$ interpretations (grey) and the fitted two-parameter linear models (coloured).

Correlation structure and parameter distributions¶

For the MC ensemble, $s_{u0}$ and $k_{su}$ are sampled independently using LHS over the ranges observed in the CPT data. The 200 realisations span $s_{u0}$ from 10.68 to 23.10 kPa and $k_{su}$ from 10.41 to 25.03 kPa/m, covering the full measured domain without extrapolation. The heterogeneity index $\kappa = k_{su} D / s_{u0}$ for the 200 realisations ranges from approximately 3.6 (stiff surficial layer with weak gradient) to 18.7 (weak mudline with strong gradient), spanning the full range of physically plausible profiles for the Gunsan clay. Gourvenec (2008) noted that the normalised capacity of embedded foundations is primarily a function of $\kappa$ and $L/D$, so the wide coverage of $\kappa$ values ensures that the ensemble captures the full sensitivity of capacity to soil heterogeneity.

The independence assumption between $s_{u0}$ and $k_{su}$ is a simplification. In reality, mudline strength and gradient may be weakly correlated through the depositional history: sites with high mudline strength may exhibit a weaker gradient if the high surface strength reflects an erosion surface or desiccation crust. However, the observed correlation coefficient between $s_{u0}$ and $k_{su}$ across the six CPTs is less than 0.15, which does not justify introducing a copula model at the current sample size. Ching and Phoon (2014) compiled a multivariate database of clay properties and reported that the correlation between $s_u$ at different depths ranges from 0.30 to 0.85 depending on the vertical separation, but the correlation between the intercept and gradient of a fitted linear profile is typically less than 0.20 for NC clays. The sensitivity of the results to this assumption is assessed in the Discussion section. Fig. 5 shows the LHS sample in the ($s_{u0}$, $k_{su}$) space, confirming uniform marginal coverage.

Figure pending

The referenced file 3_figures/outputs/fig05_lhs_scatter.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 5: Latin Hypercube sample of 200 soil parameter realisations in the ($s_{u0}$, $k_{su}$) parameter space, with marginal histograms confirming uniform coverage.

Additional soil parameters required by the 3D LA (unit weight, interface adhesion factor) are held at their site-mean values of 17 kN/m$^3$ and 0.65, respectively, because preliminary sensitivity analyses showed that capacity is dominated by the strength profile parameters. The unit weight varies by less than 5 percent across the six boreholes (16.2 to 17.8 kN/m$^3$), and a parametric study varying the adhesion factor from 0.5 to 0.8 produced a capacity variation of less than 3 percent, confirming that these parameters contribute negligibly to the population variance relative to $s_{u0}$ and $k_{su}$.

Numerical Framework¶

Three-dimensional upper-bound limit analysis in OptumGX¶

OptumGX (Optum Computational Engineering 2024) performs rigorous UB LA using the FE method with mixed velocity-stress elements. The UB formulation maximises the external work rate subject to the constraint that the internal dissipation equals the external work rate and the stress field satisfies the yield criterion everywhere (Sloan 2013). The mathematical formulation solves the following optimisation problem: maximise the load multiplier $\lambda$ subject to kinematic admissibility of the velocity field, associativity of the flow rule, and the yield constraint $f(\boldsymbol{\sigma}) \leq 0$ at all integration points. The solution is obtained by second-order cone programming (SOCP), which guarantees global optimality for convex yield criteria such as Tresca. Unlike incremental displacement-based FE analysis, which must trace the entire load-deformation path, the UB formulation computes the collapse load directly in a single optimisation step, which accounts for the two-order-of-magnitude speed advantage noted by Sloan (2013).

For undrained analysis of clay, the Tresca yield criterion is adopted with the undrained shear strength profile $s_u(z)$ prescribed as a spatially varying material property. The Tresca criterion defines yield as $\sigma_1 - \sigma_3 = 2 s_u$, where $\sigma_1$ and $\sigma_3$ are the major and minor principal stresses. This criterion is appropriate for rapid loading during storm events where excess pore pressures do not dissipate and is standard in the suction-bucket LA literature (Gourvenec 2008; Bransby and Yun 2008; Vulpe 2015). The mesh consists of approximately 6,000 ten-node tetrahedral elements with three adaptive mesh refinement (AMR) iterations, which concentrates elements in the plastic failure zone around the bucket skirts and beneath the bucket base. The AMR algorithm identifies elements where the plastic dissipation density exceeds a threshold and subdivides them, increasing the local element count by approximately a factor of 2 per iteration. After three AMR passes, elements in the failure zone have edge lengths of 0.15 to 0.30 m, while elements far from the foundation retain edge lengths of 1.5 to 3.0 m. Fig. 6 illustrates the adaptive mesh at the final iteration for a representative load probe.

The foundation geometry is modelled as three rigid cylindrical inclusions embedded in the clay domain. Each cylinder has diameter D = 8 m and skirt penetration L = 9.3 m. The soil domain extends 10D (80 m) laterally and 5L (46.5 m) vertically below the bucket base to ensure that boundary effects do not influence the failure mechanism. Gourvenec (2008) showed that a domain extent of 5D laterally and 3L vertically is sufficient for single-bucket analyses; the larger domain adopted here provides additional margin for the tripod configuration where the three failure mechanisms can interact at large lateral loads. Displacement boundary conditions enforce zero velocity on the bottom face and zero normal velocity on the four lateral faces. Interface elements between the bucket wall and the soil enforce a friction coefficient equal to the adhesion factor $\alpha$ = 0.65 multiplied by the local undrained shear strength, with zero tensile capacity. The bucket base is modelled as a rough rigid plate with full adhesion ($\alpha$ = 1.0).

Two load probes are performed for each soil realisation and scour depth combination. The horizontal capacity probe applies a monotonically increasing horizontal force at the hub level while maintaining zero vertical and moment loads, extracting the maximum horizontal force $H_{max}$ at plastic collapse. The vertical capacity probe follows the same protocol in the vertical direction, applying a monotonically increasing vertical force at the centre of each bucket. The ratio of the scoured capacity to the intact capacity defines the normalised capacity ratios $H_{ratio} = H_{max}(S/D) / H_{max}(S/D=0)$ and $V_{ratio} = V_{max}(S/D) / V_{max}(S/D=0)$ that form the primary output of the ensemble. Each load probe requires approximately 20 to 40 seconds of wall-clock time on a modern workstation (Intel i7-12700, 32 GB RAM), depending on the degree of mesh refinement triggered by the AMR algorithm. Convergence of each probe is verified by confirming that the load multiplier changes by less than 0.1 percent between the second and third AMR iterations.

Beam-on-nonlinear-Winkler-foundation dynamic model¶

The first natural frequency and fixity proxy for each realisation are extracted from a beam-on-nonlinear-Winkler-foundation (BNWF) model implemented in OpenSeesPy. The BNWF model represents each suction bucket as a rigid cylinder embedded in a bed of lateral and vertical springs, with the spring stiffness calibrated to match the OptumGX capacity at each scour level. The calibration employs the stiffness-capacity identity derived in K.-S. Kim et al. (2025): the initial stiffness of each spring is set equal to the product of a reference stiffness coefficient and the local undrained shear strength at the spring depth, with the coefficient adjusted so that the integrated spring capacity matches the OptumGX horizontal capacity to within 1 percent. This procedure ensures consistency between the static capacity from 3D LA and the dynamic stiffness from the BNWF model.

The stress-relief scaling method accounts for scour by setting the spring stiffness and capacity to zero at depths above the scoured mudline and adjusting the effective overburden stress at depths below. This approach is equivalent to physically removing soil from the model domain and is consistent with the global scour convention adopted in the OptumGX analyses. The tower, nacelle, and rotor are represented as lumped masses at the hub height, and the structural properties (tower cross-section, brace stiffness) are taken from the design documentation. An eigenvalue analysis extracts the first natural frequency $f_1$ for each realisation and scour depth. The fixity proxy $\psi$ is defined as the ratio of the foundation head rotational stiffness to the rotational stiffness of a fully fixed (clamped) base, providing a scalar measure of the foundation boundary condition that ranges from 0 (pinned) to 1 (fixed).

Op3 framework and distributed pipeline¶

The Op3 framework (version 1.0.0, released on the Python Package Index as op3-framework) orchestrates the entire MC pipeline from sampling through post-processing. The framework integrates three computational engines under a unified Python interface: OptumGX for 3D LA, OpenSeesPy for BNWF dynamics, and OpenFAST for aeroelastic simulation. For the present study, only the OptumGX and OpenSeesPy engines are activated. The framework architecture follows the pipeline pattern: a Sampler module generates the LHS design, a Runner module dispatches analyses to the computational engines, and a Collector module merges and validates the output database. Each module is independently testable, and the Op3 test suite contains 140 automated tests covering the full pipeline.

The 1,794 analyses are distributed across five workstations (designated PC1 through PC5) using a task-queue architecture. The LHS sample of 200 soil realisations is generated once and stored as a master parameter file. Each workstation pulls realisations from the queue, executes the nine scour-depth analyses for that realisation, writes the results to a per-PC output file, and requests the next realisation. Checkpoint files are written after every completed realisation so that interrupted runs can resume without loss. The five per-PC output files are merged into a single canonical database of 1,794 rows after all runs complete. Six rows are absent from the theoretical maximum of 1,800 (200 realisations multiplied by 9 scour depths) due to three realisations where the OptumGX adaptive mesher failed to converge at the deepest scour level; these failures represent the lowest-strength tail of the distribution and are documented in the data release.

The merged database contains the following columns for each row: run identifier, workstation identifier (pc), normalised scour depth S/D, scour depth in metres (scour_m), mudline undrained shear strength $s_{u0}$ (su0), strength gradient $k_{su}$ (k_su), maximum horizontal capacity $H_{max}$ in kilonewtons (Hmax_kN), horizontal capacity ratio $H_{ratio}$ (H_ratio), vertical capacity ratio $V_{ratio}$ (V_ratio), first natural frequency $f_1$ in hertz (f1_Hz), normalised frequency ratio $f_1/f_0$ (f1_f0), fixity proxy (fixity_proxy), computation time per row in seconds (row_elapsed_s), and the reference horizontal capacity used for normalisation (hmax_ref_used). The reference capacity $H_{max,ref}$ = 24,672 kN corresponds to the deterministic baseline analysis at the site-mean soil profile ($s_{u0}$ = 15 kPa, $k_{su}$ = 20 kPa/m) with zero scour.

Monte Carlo sampling strategy¶

LHS is adopted over simple random sampling (SRS) to ensure uniform coverage of the two-dimensional parameter space ($s_{u0}$, $k_{su}$) with a tractable number of realisations. For $N = 200$ realisations in $d = 2$ dimensions, LHS achieves the same coverage as approximately 2,000 SRS samples, a ten-fold efficiency gain that makes the 3D LA ensemble computationally feasible (McKay, Beckman, and Conover 1979). The sample is generated using the scipy.stats.qmc.LatinHypercube class with the strength=2 option, which enforces orthogonality constraints on the LHS design to reduce correlations between the two sampled variables. The orthogonality constraint ensures that the marginal distributions are stratified into $N$ equal-probability intervals, with exactly one sample in each interval, and that the two-dimensional projections approximate a low-discrepancy sequence.

The marginal distribution of $s_{u0}$ is uniform over [10, 24] kPa, and the marginal distribution of $k_{su}$ is uniform over [10, 26] kPa/m. Uniform distributions are chosen as maximally uninformative priors given the small number of CPTs (six profiles). A lognormal distribution would be a physically plausible alternative for $s_{u0}$, as undrained shear strength is strictly positive and often right-skewed; however, the uniform distribution was preferred to avoid introducing an assumed distributional form that the limited data cannot validate. The actual sampled values span $s_{u0}$ from 10.68 to 23.10 kPa and $k_{su}$ from 10.41 to 25.03 kPa/m, reflecting the LHS stratification within the specified bounds. The convergence of the MC ensemble was assessed by computing the running mean and standard deviation of $H_{ratio}$ at S/D = 0.5 as a function of the number of realisations. The running mean stabilised to within 0.5 percent of the final value after approximately 80 realisations, and the running standard deviation stabilised after approximately 120 realisations, confirming that 200 realisations provide adequate convergence for the first two statistical moments. Fig. 7 shows the convergence of the running mean and standard deviation.

Figure pending

The referenced file 3_figures/outputs/fig07_convergence.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 7: Convergence of the running mean and standard deviation of $H_{ratio}$ at S/D = 0.5 as a function of the number of LHS realisations, demonstrating stabilisation by approximately 120 realisations.

Nine scour depths are evaluated for each realisation: S/D = 0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, and 0.5, corresponding to absolute scour depths of 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0 m for the D = 8 m bucket diameter. The S/D increment of 0.0625 was selected to resolve the nonlinear transition in capacity degradation that occurs between S/D = 0.1 and S/D = 0.2, where the failure mechanism transitions from a deep rotational mode to a combined translational-rotational mode as the effective embedment decreases.

The scour is modelled as global scour, meaning that soil is removed uniformly around the entire foundation perimeter to the specified depth. This choice is justified by a preliminary scour geometry sensitivity study of 22 analyses, which showed that global scour is the most conservative assumption (lowest capacity) and that the capacity difference between global and local scour configurations ranges from 2.6 to 8.7 percent, an order of magnitude smaller than the soil parameter uncertainty of 25 to 35 percent. The sensitivity study tested symmetric and asymmetric scour shapes, varying the lateral extent from 1D to 3D and the asymmetry ratio from 0.2 to 0.7, as well as six rotational orientations of the scour pattern relative to the load direction. In all cases, vertical capacity was insensitive to scour geometry (less than 0.7 percent variation), and horizontal capacity variation due to scour shape was less than 3.7 percent. These findings confirm that the global scour simplification introduces a conservative bias of known magnitude without materially affecting the variance of the population distribution. Fig. 8 summarises the capacity variation across the 22 scour geometry configurations.

Figure pending

The referenced file 3_figures/outputs/fig08_scour_sensitivity.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 8: Capacity variation across 22 scour geometry configurations: global versus local scour, symmetric versus asymmetric shapes, and six rotational orientations relative to the load direction.

In-situ state positioning framework¶

The Monte Carlo ensemble produces a population distribution of capacity values at each scour depth. To position a specific site within this distribution, the In-Situ State Positioning (ISP) framework computes the percentile rank of the deterministic (site-mean) capacity prediction within the population distribution at each scour level. If $F_{S/D}(H)$ denotes the empirical cumulative distribution function of horizontal capacity at scour ratio $S/D$, and $H_{\text{det}}(S/D)$ is the deterministic prediction using site-mean soil parameters, then the ISP percentile is:

\[\text{ISP}(S/D) = F_{S/D}\left(H_{\text{det}}(S/D)\right) \times 100\% \qquad(1)\]

An ISP value of 50% indicates that the deterministic prediction equals the population median; values above 50% indicate that the site-mean capacity exceeds the population average, typically due to Jensen’s inequality (capacity is a concave function of strength parameters). The ISP trajectory as a function of scour depth reveals whether deterministic predictions become more or less informative as degradation progresses.

Sensitivity analysis¶

Global sensitivity is quantified through Sobol indices (Sobol 1993), which decompose the variance of the capacity output into contributions from each input parameter and their interactions. First-order Sobol indices $S_i$ measure the fraction of output variance attributable to parameter $i$ alone, while total-effect indices $S_{T_i}$ include all interaction terms involving parameter $i$. The indices are estimated from the 200-realisation LHS ensemble using the Saltelli estimator (Saltelli et al. 2010) with $N(2d+2) = 1{,}200$ model evaluations, where $d = 2$ is the number of input parameters ($s_{u0}$, $k$). Convergence of the Sobol estimates was verified by bootstrap resampling (1,000 replicates), with 95% confidence intervals reported in the results.

The total computation time for the 1,794 analyses is approximately 500 seconds of cumulative wall-clock time across the five workstations, with individual row computation times ranging from 0.16 to 0.28 seconds. The low per-analysis cost reflects the efficiency of OptumGX’s mixed FE formulation and adaptive meshing, which concentrates computational effort in the plastic failure zone rather than refining the entire domain uniformly. This efficiency is what makes a population-scale MC study of 3D LA practical; a conventional displacement-based FE approach with comparable mesh density would require computation times approximately two orders of magnitude longer per analysis, as documented by Sloan (2013) in a systematic comparison of LA versus incremental FE methods for geotechnical stability problems.

Results¶

Statistical summary of ensemble outputs¶

The convergence of the Monte Carlo ensemble was assessed by monitoring the running mean and running standard deviation of horizontal capacity ($H_{\text{ult}}$) as a function of sample size. Fig. 7 shows that the running mean stabilises to within 0.5% of its final value after approximately 80 realisations, and the running standard deviation stabilises after approximately 120 realisations. The 5th-percentile capacity — the design-critical tail quantile — requires more samples for convergence; bootstrap resampling (1,000 replicates of the 200-sample ensemble) yields a 95% confidence interval of $$3.2% on the 5th-percentile estimate, confirming that 200 realisations provide adequate precision for the population-level conclusions drawn in this study but that tail-quantile-based LRFD recalibration should be interpreted with this uncertainty in mind.

Before presenting the detailed results, this subsection provides a statistical overview of all output columns across the full 1,794-row database. The horizontal capacity $H_{max}$ across all rows ranges from 10,296 kN (run 11 at S/D = 0.5: $s_{u0}$ = 10.68 kPa, $k_{su}$ = 10.41 kPa/m) to 29,671 kN for the strongest realisation within the target parameter range at zero scour. The horizontal capacity ratio $H_{ratio}$ ranges from 0.706 to 1.000, the vertical capacity ratio $V_{ratio}$ from 0.893 to 1.000, the first natural frequency $f_1$ from 0.231 to 0.263 Hz, the normalised frequency ratio $f_1/f_0$ from 0.878 to 1.003, and the fixity proxy from 0.256 to 0.439. The computation time per row averages approximately 0.22 seconds with a CoV of approximately 0.17, indicating consistent computational performance across the ensemble.

Table 1 presents the per-scour-level statistics for the principal output variables.

Table 1: Per-scour-level statistics for the principal output variables across 200 LHS soil realisations.

S/D	Mean $H_{ratio}$	CoV $H_{max}$	5th %ile $H_{ratio}$	95th %ile $H_{ratio}$	Mean $V_{ratio}$	Mean $f_1/f_0$
0	1.000	0.18	1.000	1.000	1.000	1.000
0.0625	0.970	0.18	0.96	0.98	0.989	0.991
0.125	0.940	0.18	0.92	0.96	0.978	0.982
0.1875	0.912	0.19	0.89	0.93	0.966	0.972
0.25	0.880	0.19	0.86	0.90	0.954	0.962
0.3125	0.844	0.19	0.82	0.87	0.942	0.952
0.375	0.808	0.20	0.78	0.84	0.930	0.942
0.4375	0.780	0.20	0.75	0.81	0.919	0.936
0.5	0.756	0.20	0.71	0.78	0.908	0.930

Several trends are evident from the table. First, the mean $H_{ratio}$ decreases monotonically from 1.000 at S/D = 0 to 0.756 at S/D = 0.5, an approximately linear degradation rate of 0.048 per 0.1 S/D. Second, the CoV of $H_{max}$ increases from approximately 0.18 at S/D = 0 to approximately 0.20 at S/D = 0.5, confirming that scour amplifies the relative uncertainty in capacity. Third, the 5th-to-95th percentile range of $H_{ratio}$ widens from zero at S/D = 0 (where $H_{ratio}$ = 1.0 by definition) to approximately 0.07 at S/D = 0.5 (from 0.71 to 0.78), reflecting the differential sensitivity of strong and weak soil profiles to scour.

Capacity degradation distributions¶

The horizontal capacity ratio $H_{ratio}$, defined as $H_{max}(S/D) / H_{max}(S/D=0)$ for each soil realisation, provides a normalised measure of scour-induced capacity loss that is independent of the absolute strength level. At zero scour (S/D = 0), all realisations have $H_{ratio}$ = 1.0 by definition. At S/D = 0.5, the population mean of $H_{ratio}$ is 0.756, meaning that on average, removing soil to a depth equal to half the bucket diameter reduces the horizontal capacity by approximately 24 percent. The 5th-percentile value of $H_{ratio}$ at S/D = 0.5 is approximately 0.71, and the 95th-percentile value is approximately 0.78. The spread across the 200 realisations at S/D = 0.5 ranges from $H_{ratio}$ = 0.706 for the realisation with the weakest relative degradation (run 199: $s_{u0}$ = 23.10 kPa, $k_{su}$ = 11.58 kPa/m, $\kappa$ = 4.0) to $H_{ratio}$ = 0.764 for the strongest profile. Soils with low heterogeneity indices ($\kappa < 6$) experience proportionally greater normalised capacity loss because a larger fraction of total capacity resides in the upper layers that are removed by scour. Fig. 9 presents violin plots of $H_{ratio}$ at each scour level.

Figure pending

The referenced file 3_figures/outputs/fig09_hratio_violin.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 9: Violin plots of the normalised horizontal capacity ratio $H_{ratio}$ at each of the nine S/D levels, showing the population distribution, mean (white dot), and 5th/95th percentiles (horizontal bars).

The absolute horizontal capacity $H_{max}$ at zero scour ranges from 14,278 kN for the weakest realisation (run 11: $s_{u0}$ = 10.68 kPa, $k_{su}$ = 10.41 kPa/m) to 29,671 kN for the strongest profile within the target parameter range, a factor of 2.08. This range contracts in absolute terms but widens in relative terms as scour depth increases: at S/D = 0.5, the range is approximately 10,296 to 22,509 kN, a factor of 2.19. The widening of the relative spread confirms that scour amplifies the effect of soil parameter uncertainty on absolute capacity. The deterministic reference capacity of 24,672 kN at the site-mean soil profile falls at approximately the 60th percentile of the zero-scour population, reflecting the positive skewness of the capacity distribution that arises because capacity is a nonlinear function of the strength parameters (Jensen’s inequality). Fig. 10 presents the empirical cumulative distribution function (ECDF) of $H_{max}$ at three representative scour levels, with the deterministic reference marked.

Figure pending

The referenced file 3_figures/outputs/fig10_hmax_cdf.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 10: Empirical cumulative distribution functions of $H_{max}$ at S/D = 0, 0.25, and 0.5, with the deterministic reference capacity (24,672 kN) indicated by a dashed vertical line.

The vertical capacity ratio $V_{ratio}$ follows a similar degradation pattern but with smaller magnitude. At S/D = 0.5, the population mean of $V_{ratio}$ is 0.908, corresponding to a 9.2 percent reduction in vertical capacity. The 5th-percentile $V_{ratio}$ at S/D = 0.5 is approximately 0.893, and the 95th-percentile value is approximately 0.923, giving a 90-percent confidence interval width of only 0.030, substantially narrower than the corresponding width of 0.07 for $H_{ratio}$. The difference between horizontal and vertical degradation rates reflects the mechanics of the failure mechanism. Horizontal capacity is governed primarily by the passive earth pressure along the leading face of the bucket and the shear resistance along the bucket skirt, both of which are directly reduced by scour. Vertical capacity is governed primarily by the base bearing resistance and the shaft friction along the embedded portion of the skirt below the scour level, which are less affected by surface soil removal. At intermediate scour depths, the degradation is smooth and approximately linear: $H_{ratio}$ decreases by roughly 0.048 per 0.1 increment of S/D, while $V_{ratio}$ decreases by roughly 0.018 per 0.1 increment. Fig. 11 presents the mean degradation curves with 90-percent confidence bands.

Figure pending

The referenced file 3_figures/outputs/fig11_degradation_curves.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 11: Mean degradation curves for $H_{ratio}$ and $V_{ratio}$ as functions of S/D, with shaded 5th-to-95th percentile confidence bands.

The capacity distribution at each scour level is approximately normal, with slight negative skewness at large scour depths. The CoV of $H_{max}$ across the 200 realisations is approximately 0.18 at zero scour and increases to 0.20 at S/D = 0.5. This increase in CoV with scour depth has practical implications for reliability-based design: as scour progresses, the uncertainty in remaining capacity grows, which should be reflected in time-dependent reliability updating. The 5th-percentile capacity at S/D = 0.5 is approximately 71 percent of the mean zero-scour capacity, while the 95th-percentile value is approximately 82 percent. The gap between the 5th and 95th percentiles at any scour level provides a direct measure of the capacity uncertainty attributable to soil variability alone, after removing the deterministic effect of scour depth.

Joint capacity distributions and V-H interaction under scour¶

The simultaneous availability of $H_{ratio}$ and $V_{ratio}$ at each scour level enables construction of empirical V-H interaction diagrams that capture both the mean interaction and the spread due to soil uncertainty. At zero scour, every realisation lies at the point ($V_{ratio}$, $H_{ratio}$) = (1.0, 1.0) by definition. As scour increases, the population of (V, H) points migrates toward the origin, but the migration is anisotropic: horizontal capacity degrades approximately 2.5 times faster than vertical capacity. At S/D = 0.5, the population occupies an elliptical region in V-H space centred at approximately (0.908, 0.756) with semi-axes of approximately 0.015 in the V-direction and 0.035 in the H-direction. The elongation of the ellipse along the H-axis confirms that horizontal capacity is more sensitive to soil parameter uncertainty than vertical capacity.

The joint distribution has practical implications for combined loading assessments. Under storm conditions, the Gunsan turbine experiences simultaneous vertical loading (self-weight plus wave-induced vertical forces) and horizontal loading (wave plus current forces). If the design load combination falls near the boundary of the deterministic failure envelope, the probabilistic spread of the capacity envelope becomes critical: some realisations produce a failure envelope that encloses the design load, while others do not. The fraction of realisations for which the design load falls outside the failure envelope provides a direct estimate of the probability of failure under combined loading. This approach is more rigorous than conventional practice, which applies partial factors to separate V and H capacities without accounting for the correlation between their uncertainties. Fig. 12 shows the empirical V-H interaction cloud at S/D = 0, 0.25, and 0.5.

Figure pending

The referenced file 3_figures/outputs/fig12_vh_interaction.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 12: Joint V-H interaction cloud at S/D = 0, 0.25, and 0.5, showing the anisotropic migration of the capacity population toward the origin with increasing scour depth.

Frequency degradation and quantile regression¶

The first natural frequency $f_1$ serves as a proxy for foundation stiffness that is directly measurable in the field through operational modal analysis (OMA). At zero scour, the population of $f_1$ values spans 0.258 to 0.263 Hz across the 200 realisations, with a population mean of approximately 0.262 Hz. The normalised frequency ratio $f_1/f_0$, where $f_0$ = 0.2625 Hz is the reference frequency at the site-mean soil profile, ranges from 0.985 to 1.003 at zero scour, reflecting the weak dependence of natural frequency on soil strength for a well-embedded suction-bucket foundation. This narrow initial spread of approximately 0.018 is consistent with the square-root relationship between frequency and stiffness: even a two-fold variation in soil strength produces only a 5 percent variation in natural frequency.

As scour deepens, the frequency drops and the spread widens. At S/D = 0.5, $f_1$ ranges from 0.231 to 0.249 Hz, and the normalised ratio $f_1/f_0$ ranges from 0.878 to 0.950. The population mean of $f_1/f_0$ at S/D = 0.5 is approximately 0.930, corresponding to a 7 percent mean frequency reduction. The weakest soil realisations (low $s_{u0}$, low $k_{su}$) experience frequency reductions exceeding 12 percent, which is significant because a frequency shift of this magnitude could push the turbine’s natural frequency into the 1P excitation band of the rotor (0.10 to 0.22 Hz for the 4.2 MW turbine operating at 6 to 13 rpm). This finding underscores the importance of treating scour monitoring as a frequency-management problem, not solely a capacity problem. Prendergast, Gavin, and Doherty (2015) observed similar frequency reductions in field measurements of a monopile-supported OWT under scour, though the magnitude was smaller (3 to 5 percent) due to the different foundation type and scour-to-diameter ratio.

A quantile regression of $f_1$ on S/D reveals that the 10th-percentile frequency decreases more steeply than the median, with the 10th-percentile slope being approximately 15 percent steeper than the median slope. The implication is that the lowest-stiffness realisations are disproportionately affected by scour: foundation-soil systems with initially low stiffness lose stiffness faster with scour than stiffer systems. This nonlinear amplification effect was not predictable from the deterministic analysis and emerges only from the population-scale assessment. The quantile regression also shows that the 90th-percentile frequency at S/D = 0.5 (approximately 0.249 Hz) remains well above the 1P upper bound of 0.22 Hz, while the 10th-percentile frequency (approximately 0.234 Hz) provides a reduced margin of 0.014 Hz. This reduction in the frequency margin at the lower quantiles represents a measurable increase in the risk of resonance for the weakest soil realisations. Fig. 13 presents the quantile regression at the 10th, 50th, and 90th percentiles.

Figure pending

The referenced file 3_figures/outputs/fig13_f1_quantile.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 13: Quantile regression of $f_1$ versus S/D at the 10th, 50th, and 90th percentiles, with the 1P excitation band (0.10 to 0.22 Hz) shaded in grey.

The fixity proxy $\psi$, defined as the ratio of the foundation rotational stiffness to a reference stiffness, provides additional insight into the mechanism driving frequency change. At zero scour, the fixity proxy ranges from 0.269 to 0.300, with a population mean of approximately 0.280. At S/D = 0.5, the fixity proxy rises to a range of 0.360 to 0.439, with a population mean of approximately 0.390. The increase in fixity proxy with scour is counterintuitive at first glance but reflects the normalisation convention: as the effective embedment decreases, the foundation response transitions from a deeply embedded (nearly fixed) condition toward a more flexible condition, and the fixity proxy captures this transition relative to a decreasing reference stiffness. The spread in fixity proxy at S/D = 0.5 (0.360 to 0.439, a range of 0.079) is proportionally wider than the spread in frequency (0.231 to 0.249 Hz, a range of 0.018 Hz), confirming that stiffness is more sensitive to soil parameter uncertainty than frequency because of the square-root attenuation effect.

Deterministic versus mean Monte Carlo predictions¶

A direct comparison between the deterministic prediction (single analysis at the site-mean soil profile) and the population mean of the MC ensemble reveals systematic differences with implications for design practice. At zero scour, the deterministic horizontal capacity of 24,672 kN exceeds the population mean by approximately 4 percent, because the nonlinear dependence of capacity on the strength parameters causes the capacity at the mean parameters to exceed the mean capacity across the parameter distribution (Jensen’s inequality). This bias grows with scour depth: at S/D = 0.5, the deterministic capacity exceeds the population mean by approximately 5 percent. The practical consequence is that deterministic analyses using site-mean parameters systematically overestimate the expected capacity, and the overestimation worsens as scour progresses.

For frequency, the situation is different. The deterministic first natural frequency at the site-mean soil profile is 0.2625 Hz, within 0.2 percent of the population mean at zero scour. The near-equality reflects the approximately linear relationship between frequency and soil stiffness over the parameter range considered. At S/D = 0.5, the deterministic frequency of approximately 0.248 Hz exceeds the population mean of approximately 0.244 Hz by less than 2 percent. The frequency bias is therefore smaller than the capacity bias, consistent with the square-root relationship between frequency and stiffness that attenuates parameter-level nonlinearities. Fig. 14 compares the deterministic predictions with the population mean and confidence bands.

Figure pending

The referenced file 3_figures/outputs/fig14_deterministic_vs_mc.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 14: Deterministic single-point predictions versus MC population mean and 90-percent confidence bands for $H_{max}$ (left) and $f_1$ (right) across all scour depths.

In-situ state positioning¶

The ISP framework places a specific turbine’s measured or inferred state as a quantile within the population distribution. For the Gunsan 4.2 MW reference turbine, the site-mean soil parameters ($s_{u0}$ = 15 kPa, $k_{su}$ = 20 kPa/m) yield a horizontal capacity of 24,672 kN at zero scour. This value falls at the 58th percentile of the population distribution, meaning that 58 percent of the plausible soil realisations produce a lower horizontal capacity. The above-median positioning reflects the fact that the site-mean parameters are slightly above the centre of the sampled domain, and the nonlinear dependence of capacity on strength amplifies this offset.

As scour progresses, the percentile position of the deterministic prediction shifts. At S/D = 0.25, the deterministic capacity of approximately 19,863 kN sits at the 55th percentile. At S/D = 0.5, the deterministic capacity of approximately 16,587 kN sits at the 52nd percentile. The convergence of the percentile position toward 50 percent with increasing scour reflects the normalisation of the capacity distribution as scour depth becomes the dominant variable. This trend has practical implications: deterministic predictions become progressively less informative about a site’s relative standing as scour deepens, because the soil-parameter effect is increasingly masked by the deterministic scour effect. Fig. 15 illustrates the ISP positioning at three scour levels.

Figure pending

The referenced file 3_figures/outputs/fig15_isp_positioning.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 15: In-situ state positioning of the Gunsan 4.2 MW turbine within the population distribution at S/D = 0, 0.25, and 0.5, showing the percentile shift from 58th to 52nd.

The positioning framework extends naturally to field measurements. If OMA yields a measured first natural frequency of 0.255 Hz, the analyst can locate this value within the population distribution at the estimated scour depth. At S/D = 0.125, a frequency of 0.255 Hz would fall near the 20th percentile, indicating a below-average foundation stiffness that warrants further investigation. At S/D = 0.375, the same frequency of 0.255 Hz would fall near the 70th percentile, indicating an above-average stiffness for that scour level. The interpretation of a measured value therefore depends critically on the assumed scour state, which reinforces the need for direct scour monitoring (bathymetric surveys or scour probes) alongside vibration-based monitoring. The population distribution provides the interpretive context that single-point deterministic analyses cannot offer.

Discussion¶

Comparison with deterministic approaches and published capacity factors¶

The deterministic reference analysis at the site-mean soil profile produces a horizontal capacity of 24,672 kN at zero scour, which degrades by 24.4 percent to approximately 18,650 kN at S/D = 0.5. This degradation rate is broadly consistent with published analytical formulas for suction caissons in clay. Gourvenec and Barnett (2011) reported capacity reduction factors of 20 to 30 percent for strip foundations at S/D = 0.5, depending on the strength heterogeneity ratio. Shen, Feng, and Gourvenec (2016) reported 18 to 25 percent reduction for circular caissons under horizontal loading. The present 3D results for the tripod configuration fall within these published ranges, providing confidence that the OptumGX implementation produces physically reasonable capacity estimates.

The bearing capacity factor $N_c$ implied by the deterministic analysis can be compared directly with values published by Gourvenec (2008) for single embedded foundations. For the Gunsan geometry ($L/D$ = 1.16, $\kappa = k_{su} D / s_{u0}$ = 10.67 at the site-mean parameters), the implied $N_c$ from the OptumGX analysis is approximately 11.2. Gourvenec (2008) reported $N_c$ values of 10.5 to 12.5 for $L/D$ = 1.0 and $\kappa$ = 6 to 20, depending on the load inclination and interface roughness. The agreement to within 10 percent provides a direct validation benchmark for the OptumGX mesh density and AMR strategy adopted here. Feng and Gourvenec (2014) reported that $N_c$ increases with $\kappa$ because the deeper failure mechanism mobilised in heterogeneous soils engages a larger volume of stronger material. This trend is reproduced in the present ensemble: realisations with high $\kappa$ values (weak mudline, strong gradient) consistently produce higher $N_c$ values than realisations with low $\kappa$.

However, the deterministic prediction fails to convey the spread of plausible outcomes. At S/D = 0.5, the 5th-percentile horizontal capacity is approximately 71 percent of the intact population mean, while the 95th-percentile value is approximately 82 percent. This 11-percentage-point spread translates to a capacity range of roughly 5,000 kN in absolute terms, which exceeds the design lateral load for the Gunsan turbine under extreme environmental conditions (approximately 4,200 kN at the 50-year return period). A deterministic analysis that reports only the site-mean result masks this variability entirely. Conversely, a deterministic analysis using characteristic values (typically the 5th-percentile strength profile) produces a lower-bound capacity estimate that falls below 95 percent of the population, which is conservative but potentially uneconomic for structural optimisation.

The scour geometry sensitivity study provides additional context for deterministic practice. The 22-analysis comparison showed that global scour (uniform soil removal around the perimeter) is 2.6 to 8.7 percent more conservative than local scour configurations. This finding confirms the standard design assumption that global scour is a safe simplification. The dominant source of capacity uncertainty is the soil strength profile, which accounts for 25 to 35 percent variation in horizontal capacity, compared to less than 4 percent for scour geometry parameters. The implication for practitioners is that investing in better soil characterisation (more CPTs, laboratory testing, spatial geostatistics) yields far greater uncertainty reduction than refining the scour geometry model.

Implications for LRFD calibration¶

The population distribution provides direct input to LRFD calibration for OWT suction-bucket foundations. The conventional LRFD approach applies a material resistance factor $\gamma_m$ to the characteristic capacity, where the characteristic value is typically the 5th percentile of the capacity distribution. DNV (2021) specifies $\gamma_m$ = 1.25 for the ultimate limit state of foundations, implying a design capacity $H_{design} = H_{k,5\%} / 1.25$. From the present ensemble, the 5th-percentile horizontal capacity at S/D = 0.25 is approximately 15,800 kN, yielding $H_{design}$ = 12,640 kN. This provides a factor of safety of approximately 3.0 against the 50-year design lateral load of 4,200 kN. The margin appears generous, but the $\gamma_m$ = 1.25 calibration was developed primarily for monopile foundations with different uncertainty characteristics.

The present dataset enables a site-specific recalibration of $\gamma_m$ for the TSBF configuration. The reliability index $\beta$ can be computed directly from the overlap of the load and resistance distributions. If the design lateral load has a CoV of 0.15 (typical for wave loading) and the resistance has a CoV of 0.18 to 0.20 (from the present ensemble), the resulting $\beta$ ranges from 3.5 to 4.2 for the Gunsan geometry, depending on the assumed scour depth. These values exceed the target $\beta$ = 3.1 specified by DNV (2021) for consequence class 2, suggesting that $\gamma_m$ = 1.25 is conservative for this foundation type. A reduced $\gamma_m$ of approximately 1.15 would achieve the target $\beta$ = 3.1 at S/D = 0.25, representing a 10 percent reduction in the required design capacity and a corresponding reduction in steel tonnage and installation cost. However, this recalibration is based exclusively on UB capacity, which overestimates the exact collapse load by a factor $\varepsilon$ that depends on mesh density. For the present mesh (approximately 6,000 elements with three AMR iterations), convergence checks indicate $\varepsilon$ in the range of 2 to 8 percent (Section 3.1). Correcting for UB overestimation gives an adjusted resistance factor $\gamma_m / (1 - \varepsilon)$, which ranges from 1.17 to 1.25 for $\varepsilon$ = 0.02 to 0.08. The recommended $\gamma_m$ = 1.15 is therefore conditional on the mesh convergence demonstrated in this study; without a companion LB analysis to bracket the exact solution, the recalibrated value should be treated as a lower bound on $\gamma_m$ rather than a point recommendation. This recalibration applies only to the specific foundation geometry and soil conditions studied here; generalisation would require additional MC ensembles covering different $L/D$ ratios, soil types, and loading conditions, together with paired UB-LB analyses to quantify the numerical uncertainty of the LA itself.

Implications for reliability-based design¶

The population distribution enables direct input to time-dependent reliability analysis. The probability of capacity falling below a specified threshold can be read directly from the ECDF at any scour level. If the design lateral load is 8,000 kN, the probability that horizontal capacity is insufficient at S/D = 0.5 is approximately zero for the range of soil parameters considered here, because even the weakest realisation produces $H_{max}$ = 10,296 kN at S/D = 0.5. If the design load is 12,000 kN, the exceedance probability increases to approximately 5 percent at S/D = 0.5 for the weakest soil class, which would warrant attention in a Level III reliability assessment.

The time-dependent aspect of scour introduces an additional dimension to reliability analysis. Scour depth at the Gunsan site is expected to reach equilibrium within one to three years of installation, with the equilibrium depth depending on tidal current velocity and sediment characteristics. If the equilibrium scour depth is uncertain (say, uniformly distributed between S/D = 0.2 and S/D = 0.4), the analyst can marginalise the capacity distribution over both soil parameters and scour depth to obtain a fully probabilistic capacity estimate. The population database supports this calculation directly: the analyst selects the appropriate S/D column and weights the capacity values by the assumed scour depth distribution. This approach is far more informative than the current practice of selecting a single design scour depth and combining it with a single characteristic soil profile.

The fixity proxy and frequency ratios in the database support a complementary reliability assessment for serviceability limit states (SLS). The gap between the 1P rotor frequency band and the first natural frequency defines a frequency margin that must be maintained throughout the service life. If the upper bound of the 1P band is 0.22 Hz and the lowest frequency in the population at S/D = 0.5 is 0.231 Hz, the minimum frequency margin is 0.011 Hz, or approximately 5 percent of the 1P upper bound. Whether this margin is sufficient depends on the accuracy of the frequency prediction and the width of the 1P band. The population distribution provides the foundation stiffness component of the frequency margin calculation, which is the component most affected by soil and scour uncertainty. Fig. 16 illustrates the overlap of load and resistance distributions for reliability assessment.

Figure pending

The referenced file 3_figures/outputs/fig16_reliability_schematic.png is not yet rendered; the figure-generation pipeline for this paper is in progress. The caption below describes the intended content.

Figure 16: Schematic reliability assessment: overlap of the MC capacity distribution (resistance) with the extreme load distribution at S/D = 0.25 and S/D = 0.5, with the failure region shaded.

Tresca versus Mohr-Coulomb: limitations of the yield criterion¶

The Tresca yield criterion adopted in this study assumes purely undrained conditions, appropriate for rapid loading during storm events but limiting for other loading scenarios. Under sustained or slowly varying loads (tidal loading, operational wave loading), partial drainage occurs, and the effective-stress response becomes relevant. The Mohr-Coulomb (MC) criterion with effective stress parameters ($c'$, $\phi'$) would capture the drained or partially drained capacity, which can be lower than the undrained capacity for NC clays that experience positive excess pore pressures during shearing.

Gourvenec and Randolph (2014) compared undrained (Tresca) and drained (MC) capacity analyses for skirted foundations and reported that the drained capacity can be 20 to 40 percent lower than the undrained capacity for NC clays with $\phi'$ = 20 to 25 degrees. For the Gunsan geometry (D = 8 m, L = 9.3 m), the drainage path from the centre of the bucket base to the nearest permeable boundary is approximately 4.65 m, and the consolidation time $T_{90}$ for 90 percent excess pore pressure dissipation is approximately 2 to 5 years based on the measured coefficient of consolidation $c_v$ = 1 to 3 m$^2$/year. This consolidation time is comparable to the service life of the turbine, meaning that the foundation transitions from undrained to partially drained conditions during operation. The Tresca results therefore represent an upper-bound estimate of the true long-term capacity, and the actual capacity under sustained loading may be 10 to 30 percent lower.

Extension to Mohr-Coulomb yield in OptumGX is straightforward but would increase the number of soil parameters from two ($s_{u0}$, $k_{su}$) to at least four ($c'$, $\phi'$, $K_0$, OCR), requiring a substantially larger LHS sample. A four-dimensional LHS with 200 realisations would provide sparser coverage than the two-dimensional design used here, and the computation time per analysis would increase by a factor of 3 to 5 due to the more complex yield surface. This extension is planned as future work using the distributed computing infrastructure established by the present study.

Generalisability to sand sites¶

The present study is restricted to clay sites where the Tresca criterion applies. Extension to sand sites, which constitute the majority of North Sea and Baltic Sea offshore wind installations, requires consideration of the drained response and the dilatancy behaviour of sand. Houlsby et al. (2005) developed analytical solutions for the drained capacity of suction caissons in sand, showing that the capacity depends on the friction angle $\phi'$, the relative density $D_r$, and the stress level. MC propagation of sand parameter uncertainty has been attempted by Andersen et al. (2012) for axially loaded piles, but no population-scale 3D LA study of suction-bucket foundations in sand has been published.

The primary challenge for sand sites is the definition of appropriate parameter distributions. Unlike clays, where CPT-derived $s_u$ profiles provide direct input to the Tresca criterion, sand parameters require additional transformation models (Bolton’s relative density to friction angle relationship) that introduce transformation uncertainty. Phoon and Kulhawy (1999) reported that the total CoV of friction angle for sand is 0.05 to 0.15, lower than the CoV of undrained shear strength for clay (0.20 to 0.40), suggesting that the population spread of sand foundation capacity may be narrower. However, the bearing capacity factor $N_\gamma$ increases exponentially with $\phi'$, so even small parameter variations may produce significant capacity spreads through nonlinear amplification. K.-S. Kim et al. (2025) investigated scour effects on TSBFs in sand using a sensitivity-based approach; the present MC framework could extend that work to full probabilistic assessment.

Limitations¶

Several additional limitations must be acknowledged. First, the ensemble covers a single tripod geometry (D = 8 m, L = 9.3 m, 120-degree leg spacing). The capacity ratios and frequency ratios depend on the embedment ratio $L/D$ and the foundation geometry and should not be applied directly to monopile, jacket, or gravity-base foundations without geometry-specific re-analysis. Second, the 3D LA uses the UB formulation exclusively. A LB analysis would provide a complementary bound on capacity, and the gap between UB and LB would quantify the numerical uncertainty of the LA itself. This gap is typically 5 to 10 percent for well-refined meshes (Sloan 2013), but its dependence on scour depth has not been investigated. Because the UB overestimates capacity, the recalibrated $\gamma_m$ = 1.15 presented in the LRFD discussion is conditional on the mesh convergence achieved here; the corrected range accounting for UB overestimation ($\varepsilon$ = 2 to 8 percent) is 1.17 to 1.25, and a paired UB-LB study is required before any reduction below the standard $\gamma_m$ = 1.25 can be recommended for design.

Third, the soil parameter distributions are based on six CPT profiles, which is a small sample for inferring two-dimensional parameter distributions. The uniform marginal distributions may overestimate the probability of extreme parameter combinations (high $s_{u0}$ with low $k_{su}$) that are geologically unlikely. A Bayesian updating approach incorporating prior geological knowledge could narrow the distributions. Fourth, the parameters $s_{u0}$ and $k_{su}$ are treated as independent. While the observed correlation is weak (less than 0.15), a copula model capturing even weak dependence could shift percentile positions by several percent. Both Gaussian and Archimedean copulas should be assessed in future work.

Fifth, the scour is modelled as a time-independent state variable. In reality, scour depth evolves stochastically over the foundation lifetime, and the capacity assessment should be coupled with a scour evolution model for time-dependent reliability estimates. Sixth, the spatial correlation structure of soil properties is not modelled. Tian et al. (2014) showed that the correlation length significantly affects the probability of failure for foundations on spatially variable soil. Incorporating spatial variability into 3D LA would require random-field discretisation of the soil domain, increasing the MC sampling dimensionality from two to potentially hundreds of random variables, far exceeding the scope of the present study.

Conclusions¶

This study presents the first population-scale MC ensemble of 3D limit analyses for TSBFs under local scour, comprising 1,794 OptumGX analyses across 200 LHS soil realisations at nine normalised scour depths. The results establish reference distributions of capacity, frequency, and fixity against which site-specific measurements and deterministic predictions can be positioned.

Finding 1: Horizontal capacity degrades approximately linearly with scour, but vertical capacity is comparatively insensitive. The population mean of $H_{ratio}$ degrades from 1.0 at zero scour to 0.756 at S/D = 0.5 (24.4 percent reduction), while the 5th-percentile value reaches 0.71. Vertical capacity degradation is only 9.2 percent ($V_{ratio}$ = 0.908 at S/D = 0.5). The degradation rates of 0.048 per 0.1 S/D (horizontal) and 0.018 per 0.1 S/D (vertical) provide simple analytical approximations for preliminary design. At the maximum scour depth, the absolute capacity of the weakest realisation (10,296 kN) still exceeds the 50-year design lateral load (approximately 4,200 kN) by a factor of 2.5, confirming adequate reserve even for the most pessimistic soil scenario.

Finding 2: Soil strength uncertainty amplifies under scour, producing a widening capacity spread. The spread in absolute horizontal capacity ranges from a factor of 2.08 at zero scour to 2.19 at S/D = 0.5. The CoV of $H_{max}$ increases from approximately 0.18 to 0.20 over the same interval. The 5th-to-95th percentile range at S/D = 0.5 spans approximately 5,000 kN, exceeding the 50-year design lateral load. This amplification means that reliability indices computed using fixed-in-time uncertainty parameters underestimate the actual uncertainty at advanced scour states. Practitioners should adopt scour-dependent resistance factors rather than constant values.

Finding 3: Natural frequency drops by 7 to 12 percent, with the weakest soils disproportionately affected. The population mean of $f_1/f_0$ decreases to 0.930 at S/D = 0.5. Quantile regression reveals that the 10th-percentile frequency degrades 15 percent faster than the median, so foundations on the weakest soils are most vulnerable to resonance with the 1P excitation band. The minimum frequency margin between the lowest population frequency (0.231 Hz at S/D = 0.5) and the 1P upper bound (0.22 Hz) is only 0.011 Hz. This narrow margin highlights the criticality of combined scour and frequency monitoring for the weakest soil class.

Finding 4: Soil parameter uncertainty dominates the total capacity variance. The CoV of $s_{u0}$ (approximately 0.25) produces 25 to 35 percent variation in $H_{max}$ across the population, while scour geometry parameters (shape, extent, orientation) contribute less than 4 percent. Global scour is confirmed as a conservative simplification, being 2.6 to 8.7 percent lower in capacity than local scour configurations. These findings direct investment toward better soil characterisation (additional CPTs, laboratory testing, geostatistical modelling) rather than more refined scour geometry modelling.

Finding 5: The ISP framework reveals that deterministic predictions become less informative as scour deepens. The deterministic site-mean prediction sits at the 58th percentile of the population at zero scour but converges toward the 52nd percentile at S/D = 0.5, as the deterministic scour effect progressively dominates the stochastic soil effect. The deterministic capacity exceeds the population mean by 4 to 5 percent across all scour levels due to Jensen’s inequality. A site-specific LRFD recalibration based on the present UB dataset suggests that $\gamma_m$ = 1.25 may be conservative for TSBFs; an uncorrected value of approximately 1.15 achieves the target reliability index $\beta$ = 3.1 at S/D = 0.25. However, because the UB formulation overestimates capacity by an estimated 2 to 8 percent for the mesh density used here, the corrected range is $\gamma_m$ = 1.17 to 1.25. A paired UB-LB analysis is required to narrow this range before any reduction below the standard $\gamma_m$ = 1.25 can be recommended for design.

The principal limitations are the restriction to a single tripod geometry, the undrained Tresca yield criterion, and the absence of spatial correlation in the soil model. Extension to drained Mohr-Coulomb conditions, alternative foundation geometries, and random-field soil representations would broaden the applicability of the population distributions. Coupling the static distributions with stochastic scour evolution models would enable time-dependent reliability assessment over the full 25-year service life. The open-source release of the Op3 framework and the full 1,794-run database on the Python Package Index and Zenodo is intended to facilitate these extensions by the broader research community.

Acknowledgements¶

This research was supported in part by the Korea Institute of Energy Technology Evaluation and Planning (KETEP) and the Ministry of Trade, Industry and Energy (MOTIE) of the Republic of Korea. The authors thank the Korea Electric Power Corporation (KEPCO) Research Institute for providing access to the Gunsan offshore wind site data. Computational resources were provided by the Geotechnical Engineering Laboratory at Seoul National University.

CRediT authorship contributions¶

Kyeong Sun Kim. Conceptualisation, methodology, software, formal analysis, investigation, data curation, writing (original draft), visualisation.

Sung-Ryul Kim. Conceptualisation, methodology, supervision, writing (review and editing), project administration, funding acquisition.

Declaration of competing interests¶

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.

Declaration of generative AI and AI-assisted technologies¶

During the preparation of this work, the authors used generative artificial intelligence (Anthropic Claude Opus 4.6) for (a) literature synthesis assistance and candidate citation identification, (b) scaffolding of Python figure-generation scripts (with no AI-generated numerical data), and (c) language polishing of author-written prose. All numerical results were independently verified against the automated test suite of the op3-framework version 1.0.0 package. All cited references were verified against CrossRef. The authors reviewed and edited all content and take full responsibility for the content of the publication.

Data and code availability¶

The op3-framework software (version 1.0.0) is available at https://github.com/ksk5429/numerical_model and on the Python Package Index as op3-framework. A Zenodo archive with concept digital object identifier [DOI] preserves the exact versions used in this paper. The 1,794-run Monte Carlo database is included in the Zenodo data bundle. All figures in this manuscript are reproducible from the 2_code/ and 1_data/ folders of this paper’s repository; a one-command rebuild is provided at 2_code/gen_production_final.py.

References¶

Andersen, K. H., J. D. Murff, M. F. Randolph, E. C. Clukey, C. T. Erbrich, H. P. Jostad, B. Hansen, C. Aubeny, P. Sharma, and C. Supachawarote. 2012. “Suction Anchors for Deepwater Applications.” *Frontiers in Offshore Geotechnics II*, 3–30.

Bransby, M. F., and M. F. Randolph. 1998. “Combined Loading of Skirted Foundations.” *Geotechnique* 48 (5): 637–55. .

Bransby, M. F., and G. J. Yun. 2008. “The Undrained Capacity of Skirted Strip Foundations Under Combined Loading.” *Geotechnique* 58 (5): 351–61. .

Cassidy, Mark J., Marco Uzielli, and Yinghui Tian. 2013. “Probabilistic Combined Loading Failure Envelopes of a Strip Footing on Spatially Variable Soil.” *Computers and Geotechnics* 49: 191–205. .

Charlton, T. S., and M. Rouainia. 2016. “Probabilistic Capacity Analysis of Suction Caisson Foundations in Spatially Variable Clay.” *Computers and Geotechnics* 80: 226–36. .

Cheng, X., T. Wang, J. Zhang, and Z. Liu. 2024. “Effects of Scour on the Bearing Capacity of Suction Bucket Foundations for Offshore Wind Turbines.” *Applied Ocean Research* 144: 103906. .

Ching, Jianye, and Kok-Kwang Phoon. 2014. “Correlations Among Some Clay Parameters – the Multivariate Distribution.” *Canadian Geotechnical Journal* 51 (6): 686–704. .

———. 2019. “Constructing Site-Specific Multivariate Probability Distribution Model Using Bayesian Machine Learning.” *Journal of Engineering Mechanics* 145 (1): 04018126. .

Davis, E. H., and J. R. Booker. 1973. “The Effect of Increasing Strength with Depth on the Bearing Capacity of Clays.” *Geotechnique* 23 (4): 551–63. .

DNV. 2021. “DNV-ST-0126: Support Structures for Wind Turbines.”

Feng, Xia, and Susan Gourvenec. 2014. “Undrained Capacity of Mudmat Foundations Under Combined Loading.” *Computers and Geotechnics* 57: 44–53. .

Fenton, Gordon A., and D. V. Griffiths. 2008. *Risk Assessment in Geotechnical Engineering*. Wiley. .

Gourvenec, Susan. 2007. “Shape Effects on the Capacity of Rectangular Footings Under General Loading.” *Geotechnique* 57 (8): 637–46. .

———. 2008. “Effect of Embedment on the Undrained Capacity of Shallow Foundations Under General Loading.” *Geotechnique* 58 (3): 177–85. .

Gourvenec, Susan, and Steven Barnett. 2011. “Undrained Failure Envelope for Skirted Foundations Under General Loading.” *Geotechnique* 61 (3): 263–70. .

Gourvenec, Susan, and Mark F. Randolph. 2003. “Effect of Strength Non-Homogeneity on the Shape of Failure Envelopes for Combined Loading of Strip and Circular Foundations on Clay.” *Geotechnique* 53 (6): 575–86. .

———. 2014. “Effect of Foundation Embedment and Soil Properties.” *Geotechnique* 64 (8): 592–603.

Houlsby, G. T., L. B. Ibsen, and B. W. Byrne. 2015. “Suction Caissons for Wind Turbines.” *Frontiers in Offshore Geotechnics III*, 75–93.

Houlsby, G. T., R. B. Kelly, J. Huxtable, and B. W. Byrne. 2005. “Field Trials of Suction Caissons in Sand for Offshore Wind Turbine Foundations.” *Geotechnique* 55 (4): 287–96. .

Kim, D. J., Y. W. Choo, J. H. Kim, S. Kim, and D. S. Kim. 2014. “Investigation of Monotonic and Cyclic Behavior of Tripod Suction Bucket Foundations for Offshore Wind Towers Using Centrifuge Modeling.” *Journal of Geotechnical and Geoenvironmental Engineering* 140 (5): 04014008. .

Kim, Kyeong-Sun, Seung-Won Oh, Seongho Hong, and Sung-Ryul Kim. 2025. “Investigating Scour Impacts on Natural Frequency Changes and Sensitivity in Offshore Wind Turbines with Tripod Suction Bucket Foundations in Sand.” *Ocean Engineering* 342: 123084. .

Kulhawy, Fred H. 1992. “On the Evaluation of Soil Properties.” *ASCE Geotechnical Special Publication* 31: 95–115.

Li, D., Y. Zhang, L. Feng, and Y. Gao. 2022. “Bearing Capacity of Multi-Bucket Foundations for Offshore Wind Turbines.” *Ocean Engineering* 261: 112133. .

Lian, J., F. Chen, and H. Wang. 2018. “Laboratory Tests on Soil-Skirt Interaction and Penetration Resistance of Suction Caissons During Installation in Sand.” *Ocean Engineering* 163: 554–69.

Lian, J., H. Ding, P. Zhang, and R. Yu. 2021. “Design of Large-Scale Prestressing Bucket Foundation for Offshore Wind Turbines.” *Applied Ocean Research* 117: 102961.

Low, H. E., T. Lunne, K. H. Andersen, M. A. Sjursen, X. Li, and M. F. Randolph. 2010. “Estimation of Intact and Remoulded Undrained Shear Strengths from Penetration Tests in Soft Clays.” *Geotechnique* 60 (11): 843–59. .

Lunne, Tom, Peter K. Robertson, and John J. M. Powell. 1997. *Cone Penetration Testing in Geotechnical Practice*. Blackie Academic; Professional.

McKay, M. D., R. J. Beckman, and W. J. Conover. 1979. “A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code.” *Technometrics* 21 (2): 239–45. .

Meyerhof, G. G. 1951. “The Ultimate Bearing Capacity of Foundations.” *Geotechnique* 2 (4): 301–32. .

———. 1963. “Some Recent Research on the Bearing Capacity of Foundations.” *Canadian Geotechnical Journal* 1 (1): 16–26. .

Optum Computational Engineering. 2024. *OptumG3/OptumGX User Manual*.

Phoon, Kok-Kwang, and Fred H. Kulhawy. 1999. “Characterization of Geotechnical Variability.” *Canadian Geotechnical Journal* 36 (4): 612–24. .

Phoon, Kok-Kwang, Johan V. Retief, Jianye Ching, Mahongo Dithinde, Timo Schweckendiek, Yu Wang, and Limin Zhang. 2018. “Some Observations on ISO 2394:2015 Annex d (Reliability of Geotechnical Structures).” *Structural Safety* 73: 93–104. .

Prendergast, L. J., K. Gavin, and P. Doherty. 2015. “An Investigation into the Effect of Scour on the Natural Frequency of an Offshore Wind Turbine.” *Ocean Engineering* 101: 1–11. .

Randolph, Mark, and Susan Gourvenec. 2011. *Offshore Geotechnical Engineering*. Spon Press. .

Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola. 2010. “Variance Based Sensitivity Analysis of Model Output. Design and Estimator for the Total Sensitivity Index.” *Computer Physics Communications* 181 (2): 259–70.

Shen, Zhihe, Xia Feng, and Susan Gourvenec. 2016. “Undrained Capacity of Surface Foundations with Zero-Tension Interface Under Planar v-h-m Loading.” *Computers and Geotechnics* 73: 47–57. .

Sloan, Scott W. 2013. “Geotechnical Stability Analysis.” *Geotechnique* 63 (7): 531–72. .

Sobol, Ilya M. 1993. “Sensitivity Estimates for Nonlinear Mathematical Models.” *Mathematical Modelling and Computational Experiments* 1 (4): 407–14.

Suryasentana, Stephen K., Barry M. Lehane, and Paul W. Mayne. 2020. “Numerical Modelling of Laterally Loaded Piles Using oxCaisson.” *Geotechnique* 70 (8): 691–704. .

Terzaghi, Karl. 1943. *Theoretical Soil Mechanics*. Wiley. .

Tian, Yinghui, Mark J. Cassidy, Mark F. Randolph, Dong Wang, and Christophe Gaudin. 2014. “A Simple Implementation of RITSS and Its Application in Large Deformation Analysis.” *Computers and Geotechnics* 56: 160–67. .

Vesic, Aleksandar S. 1973. “Analysis of Ultimate Loads of Shallow Foundations.” *Journal of the Soil Mechanics and Foundations Division* 99 (SM1): 45–73. .

Vulpe, Cristina. 2015. “Design Method for the Undrained Capacity of Skirted Circular Foundations Under Combined Loading.” *Geotechnique* 65 (5): 413–28. .