Paper deep dive
Human-guided physics-constrained AI agents construct an auditable model of soil-plug evolution
Jie Shi, Yimin Lu, Zhongkun Ouyang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 9/27/2026, 2:36:05 AM
Summary
This paper introduces a human-in-the-loop, physics-constrained multi-agent AI workflow for constructing auditable engineering models. Applied to soil-plug evolution during suction-caisson installation, the system generated and audited six candidate formulations (V1-V6) in 2.9 hours. The selected model (V6), which integrates geometric displacement, seepage-driven void-ratio evolution, and near-wall dilation, significantly reduced prediction errors compared to baselines. The workflow demonstrates that AI agents can automate the theory-to-solver chain while maintaining physical admissibility and enabling independent auditing, extending multi-agent AI beyond simple task automation.
Entities (8)
Relation Signals (5)
Human-guided physics-constrained AI agents → appliedto → Soil-plug evolution
confidence 95% · Applied to soil-plug evolution during suction-caisson installation, the workflow generated and audited 6 formulations in 2.9 h of agent execution once physical knowledge and inputs were prepared.
V6 → incorporates → Seepage-driven soil void-ratio evolution
confidence 94% · Among these formulations, adding seepage-driven soil void-ratio evolution to the geometric baseline reduced mean absolute final-heave error from 58.4% to 9.0% across 14 profiles; the selected model further incorporated near-wall dilation...
V6 → incorporates → Near-wall dilation
confidence 93% · the selected model further incorporated near-wall dilation and achieved mean absolute percentage errors of 12.4% across 9 final-state cases and 4.2% at the endpoints of 5 process histories.
V6 → outperforms → Geometric baseline
confidence 92% · adding seepage-driven soil void-ratio evolution to the geometric baseline reduced mean absolute final-heave error from 58.4% to 9.0% across 14 profiles
Human-guided physics-constrained AI agents → enables → Auditable model construction
confidence 90% · Overall, this work extends multi-agent AI beyond task automation toward human-governed engineering solvers.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Engineering predictions require physical mechanisms to be translated consistently into equations, discretization, code, and validation, yet errors can propagate despite local checks. Artificial-intelligence (AI) agents automate scientific tasks, but coordinating and independently auditing the theory-to-solver process under physical constraints and human oversight remains unresolved. We introduce a human-in-the-loop, physics-constrained multi-agent workflow where human experts define admissible physics and modeling boundaries, while agents retrieve evidence, derive equations, implement solvers, and audit the theory-to-code chain. Applied to soil-plug evolution during suction-caisson installation, the workflow generated and audited 6 formulations in 2.9 h of agent execution once physical knowledge and inputs were prepared. Among these formulations, adding seepage-driven soil void-ratio evolution to the geometric baseline reduced mean absolute final-heave error from 58.4% to 9.0% across 14 profiles; the selected model further incorporated near-wall dilation and achieved mean absolute percentage errors of 12.4% across 9 final-state cases and 4.2% at the endpoints of 5 process histories. Beyond predictive performance, blinded replay recovered all 9 target problems, while an independent audit uncovered 5 implementation problems after 36 predefined checks had passed. Overall, this work extends multi-agent AI beyond task automation toward human-governed engineering solvers.
Tags
Links
- Source: https://arxiv.org/abs/2609.23360v1
- Canonical: https://arxiv.org/abs/2609.23360v1
Trouble viewing inline? Open PDF directly →
Full Text
230,052 characters extracted from source content.
Expand or collapse full text
Human-guided physics-constrained AI agents construct an auditable model of soil-plug evolution Jie Shi Address: Institute of Ocean Engineering, Shenzhen International Graduate School, Tsinghua University, Shenzhen 518055, China Yimin Lu Email: Yimin.Lu@ttu.edu Corresponding author: Corresponding author Address: Department of Civil, Environmental, and Construction Engineering, Texas Tech University, USA Zhongkun Ouyang Email: ouyangzk@sz.tsinghua.edu.cn Corresponding author: Corresponding author Address: Institute of Ocean Engineering, Shenzhen International Graduate School, Tsinghua University, Shenzhen 518055, China Abstract Engineering predictions require physical mechanisms to be translated consistently into equations, discretization, code, and validation, yet errors can propagate despite local checks. Artificial-intelligence (AI) agents automate scientific tasks, but coordinating and independently auditing the theory-to-solver process under physical constraints and human oversight remains unresolved. We introduce a human-in-the-loop, physics-constrained multi-agent workflow where human experts define admissible physics and modeling boundaries, while agents retrieve evidence, derive equations, implement solvers, and audit the theory-to-code chain. Applied to soil-plug evolution during suction-caisson installation, the workflow generated and audited 6 formulations in 2.9 h of agent execution once physical knowledge and inputs were prepared. Among these formulations, adding seepage-driven soil void-ratio evolution to the geometric baseline reduced mean absolute final-heave error from 58.4% to 9.0% across 14 profiles; the selected model further incorporated near-wall dilation and achieved mean absolute percentage errors of 12.4% across 9 final-state cases and 4.2% at the endpoints of 5 process histories. Beyond predictive performance, blinded replay recovered all 9 target problems, while an independent audit uncovered 5 implementation problems after 36 predefined checks had passed. Overall, this work extends multi-agent AI beyond task automation toward human-governed engineering solvers. Introduction Engineering models become actionable when physical understanding can be translated into executable and validated predictions. This need is especially acute for offshore renewable-energy infrastructure, where rapid deployment must be matched by reliable foundation design under uncertain construction and loading conditions[1, 2, 3, 4, 5]. Current practice often requires a compromise between case-specific, resource-intensive experiments or high-fidelity simulations and more efficient, generalizable analytical models [5, 6, 7]. Developing such models is not a single reasoning task, but a sequence of interdependent decisions spanning problem formulation, mechanism selection, governing-equation derivation, constitutive closure, discretization, code implementation, verification, and benchmark evaluation (Figure 1(a)) [8, 9, 10]. Each stage draws on distinct expertise, and subtle inconsistencies in assumptions, notation, boundary conditions, or state updates can propagate through the workflow and remain concealed until validation. This long and fragile theory-to-solver process is not unique to offshore engineering; it also arises in physics-intensive problems such as additive-manufacturing process design, wearable-sensor development, and advanced battery design [11, 12, 13]. Suction caissons make this general difficulty concrete. These foundations use a pressure difference across the caisson lid to drive penetration into the seabed, offering installation advantages in energy demand, noise, and soil disturbance [5, 6]. However, during installation in sand, soil entering the caisson can rise above the surrounding seabed and form a hidden internal plug. Its continued growth may impede penetration, leave the caisson short of its design depth, and introduce uncertainty into foundation performance [5, 14, 15, 16, 17]. Predicting this evolution requires the simultaneous treatment of a moving boundary, suction-induced seepage, changes in effective stress and soil state, and shear and dilatancy along the soil–caisson interface (Figure 1(b)) [18, 19]. Since these processes occur inside a buried and evolving domain, experimental observations are incomplete, high-fidelity simulations are often case-specific, and reduced-order theory is difficult to formulate and implement consistently. These challenges raise the following question: can we establish an efficient, generalizable, physics-constrained, and auditable theory-to-solver framework by leveraging existing knowledge and emerging artificial-intelligence (AI) techniques? Figure 1: The theory-to-solver bottleneck and its demonstration problem. (a) Model development in physics-intensive engineering traverses interdependent stages from problem definition and mechanism selection to governing equations, constitutive closure, discretization, code implementation, verification, and benchmark evaluation. Representative applications illustrate the associated expertise, time, and error-traceability challenges. (b) Soil-plug evolution during suction-caisson installation in sand provides the demonstration problem. The hidden response couples large deformation and a moving boundary, suction-induced seepage, effective-stress and void-ratio evolution, and soil–caisson interface behavior. Large language models (LLMs) and LLM-driven multi-agent systems offer new capabilities for coordinating such workflows through natural-language reasoning, tool use, information retrieval, and conversion of modeling logic into executable operations [20, 21, 22]. Retrieval-augmented generation (RAG) can ground model outputs in external sources [23], while skill-routed retrieval can direct the search to task-specific procedures or knowledge modules [24]. However, retrieval relevance is not equivalent to physical admissibility. Before a relation can enter a mechanics model, its provenance, variables, units, assumptions, boundary conditions, and range of validity must be explicit, and it must remain compatible with the other model components. Even a relation drawn from a credible source can yield an invalid solver if it is applied outside its admissible conditions or connected incorrectly to surrounding equations. Theory implementation and code generation pose a related but distinct challenge. Recent studies have demonstrated the potential of LLMs to translate mathematical formulations into scientific-computing code, including partial-differential-equation (PDE) solvers and geotechnical finite-element models (FEMs) [25, 26]. Multi-agent architectures further decompose scientific workflows by assigning derivation, implementation, and verification to specialized roles [27, 28]. However, this decomposition also introduces handoffs among physical assumptions, continuous equations, discrete updates, and code states. Errors at these interfaces can produce software that is syntactically correct and may even pass local tests while remaining inconsistent with the intended physics or obscuring the path from source evidence to prediction [29, 30, 31]. Therefore, the central methodological challenge is how knowledge retrieval, physical reasoning, numerical implementation, and validation can be coordinated under a common set of constraints and made auditable across the complete theory-to-solver chain. In this study, we introduce a human-in-the-loop, physics-constrained multi-agent workflow for constructing, implementing, and validating engineering models across the complete theory-to-solver sequence. Human experts define the modeling objectives, admissible physical knowledge, parameter-source rules, and operational boundaries, while specialized agents collaborate on knowledge retrieval, mechanism formulation, equation derivation, numerical implementation, and independent audit. We demonstrate the workflow by developing a soil-plug evolution model for suction-caisson installation that unifies geometric displacement, seepage-induced soil-state evolution, and soil–caisson interfacial dilatancy through solid-volume conservation. The model-development process is examined through retained workflow records and role-separated audits, while the resulting solver is evaluated using staged benchmark evidence that covers formula reduction and volume closure, final state prediction, and continuous process prediction. The results show that the workflow can produce an accurate and traceable engineering solver while identifying cross-layer inconsistencies that conventional implementation checks may overlook, providing a general framework for auditable multi-agent modeling in physics-intensive engineering. Results Workflow-guided model evolution yields the selected soil-plug formulation We applied the physics-constrained multi-agent workflow to develop 6 candidate formulations (V1–V6) for soil-plug evolution during suction-caisson installation. Each version arose from a diagnostic question prompted by the performance of the preceding formulation and was evaluated using the same historical benchmark cases and parameter-source rules. Screening based on final-heave and process final-point errors identified V6 as the selected formulation since it retained a low final-heave error while producing the lowest process final-point error (Table 1 and Fig. 2). Model development proceeded through three diagnostic cycles. Table 1: Workflow-guided development and historical screening of candidate soil-plug formulations. Version Diagnostic question Workflow action Screening outcome V1 Can geometric displacement caused by skirt penetration alone explain plug heave? Established a geometric baseline from solid-volume conservation and suction-caisson design relations [5]. Produced the largest final-heave and process final-point errors among the 6 versions. V2 How does suction-driven seepage change effective stress and void ratio? Added source-supported relations linking seepage, effective stress, and void ratio to update the soil state [32, 33, 34]. Both errors decreased only slightly relative to V1. V3 How should lateral effective stress and soil rebound evolve during confined unloading? Updated lateral and mean effective stresses using relations for earth pressure at rest, stress history, and granular-soil elasticity [35, 36, 37]. Both errors decreased only modestly relative to V2. V4 Does stress-path evolution produce additional void-ratio change in dense sand? Introduced a stress-path-dependent void-ratio update based on state-parameter and state-dependent dilatancy relations [34, 38, 39]. Produced the first substantial reduction in both screening errors relative to V3. V5 Can a more compact void-ratio update improve both final-state and process-level predictions? Tested an alternative void-ratio update under the same parameter-source rules. Reduced the process final-point error but increased the final-heave error. V6 Can the admitted mechanisms be assembled to keep both final-state and process errors low? Combined the audited seepage and void-ratio relations with near-wall soil dilation and consistent final-state updating [40, 15, 6]. Retained a low final-heave error and achieved the lowest process final-point error; V6 was selected. The first cycle expanded the model beyond geometric displacement to account for suction-induced changes in soil state. V1 used solid-volume conservation to calculate the geometric heave caused by skirt penetration but did not account for changes in soil void state under suction. The marked discrepancy from the measured plug heave indicated that geometric displacement alone could not explain the observed plug heave response. Therefore, the workflow focused the next diagnostic question on how suction-induced seepage changes the effective stress and void ratio of the soil inside the caisson. To address this question, the knowledge engine retrieved source-supported relations linking suction, seepage, effective stress and void ratio [32, 33, 34], and assembled them into V2. V2 altered the predicted plug heave only marginally. An audit of the update chain found that the vertical effective stress responded to hydraulic loading, while the lateral effective stress remained tied to a fixed stress ratio. Since the soil column inside the caisson is laterally constrained, the workflow next examined the evolution of lateral effective stress during vertical unloading. Relations describing earth pressure at rest, stress history, and the elastic response of granular soils were incorporated into V3, which updated the lateral and mean effective stresses along a confined-unloading path [35, 36, 37]. This revision reduced the final-heave error only modestly. Figure 2: Historical screening of candidate soil-plug formulations. All 6 formulations were evaluated using the same benchmark cases and parameter-source rules. (a) Final-heave mean absolute percentage error (MAPE) across 9 final-state cases. (b) Process final-point MAPE across 5 measured heave–depth histories. (c) Joint evolution of the two errors across successive formulations; the inset enlarges the differences among V4–V6. Lower values indicate closer agreement with the benchmark observations. Final-heave mean absolute percentage error (MAPE) remained above 70% for V1–V3 (Figure 2), indicating that the represented volumetric response remained insufficient to reproduce the measured plug heave. The audit localized the major deficiency to the conversion of stress change into void-ratio change. Then, the second diagnostic cycle examined how the distance from the critical state bounds the available volumetric change in dense sand and how movement along the stress path mobilizes that change. The knowledge engine combined the resulting state-parameter and state-dependent dilatancy relations to construct the stress-path-dependent void-ratio update in V4. This revision produced the first substantial improvement in both screening metrics: final-heave and process final-point MAPE fell to 12.3% and 10.3%, respectively (Figure 2). The third diagnostic cycle tested whether final-state and process-level performance could be improved simultaneously. Under the same parameter-source rules, V5 tested a more compact void-ratio update. Its process final-point MAPE decreased to 6.6%, but its final-heave MAPE increased to 15.0%, revealing a trade-off between the two evaluation targets. This trade-off motivated V6, which combined the audited seepage and void-ratio relations with near-wall soil dilation and consistent final-state updating. V6 achieved a final-heave MAPE of 12.4% and reduced the process final-point MAPE to 4.2%. Its final-heave performance remained comparable to that of V4, while its process final-point error was the lowest among the 6 formulations. This balance supported the selection of V6 for subsequent analysis. After the physical knowledge and model inputs had been prepared, development of V1–V6 required 171.5 min (2.9 h) of total agent execution time (Table 2). Of this total, derivation required 42.5 min, implementation required 82.5 min, and independent audit required 46.5 min, with implementation accounting for the largest share. Execution time for an individual formulation ranged from 19.3 min for V1 to 37.5 min for V4. Across all 6 formulations, the recorded agent used 73.07 million tokens, corresponding to an estimated API-equivalent cost of US$42.23 under the GPT-5.4 standard rates available in June 2026. These measurements cover the agent execution required to translate staged modeling questions into implemented and independently audited solver versions; they do not include the preceding preparation of physical knowledge and model inputs or the time required for human review. Table 2: Agent execution time, token use, and estimated API-equivalent cost for the V1–V6 model-development. Version Derivation (min) Implementation (min) Independent audit (min) Total (min) Tokens (million) API equivalent cost (US$) V1 4.1 8.5 6.7 19.3 2.96 2.03 V2 10.8 17.6 7.4 35.7 5.14 3.82 V3 6.0 8.6 7.3 21.9 8.79 5.57 V4 5.8 21.2 10.5 37.5 15.44 8.70 V5 6.3 7.5 7.4 21.2 18.01 9.72 V6 9.4 19.1 7.3 35.8 22.72 12.39 Total 42.5 82.5 46.5 171.5 73.07 42.23 Execution time is the sum of the derivation, implementation, and independent-audit role durations recorded after the physical knowledge and model inputs had been prepared. Token totals comprise input and output tokens across the corresponding recorded agent turns; cached input is counted once within input, and reasoning tokens are included within output. API-equivalent costs use the GPT-5.4 standard rates available in June 2026: US$2.50, US$0.25, and US$15.00 per million uncached-input, cached-input, and output tokens, respectively. Model selection remains stable across most study exclusions After screening against the full literature pool selected V6, we performed a five-fold retrospective leave-one-study-out analysis to test whether this selection was dominated by an individual source study or by the weighting assigned to endpoint and process-curve errors. The 6 candidate formulations were fixed before this analysis. In each fold, one study was excluded, V1–V6 were rescored using the remaining four studies, and the selected version was evaluated against the held-out study. The analysis was repeated using endpoint-to-curve error weights of 75:25, 50:50, and 25:75. Figure 3(a) presents the candidate scores under the primary 75:25 weighting, Figure 3(b) compares the selections across all three weightings, and Figure 3(c) reports the held-out errors of the selected versions. Figure 3: Retrospective leave-one-study-out assessment of candidate-selection stability. The 6 candidate formulations were fixed before the analysis. In each fold, one study was excluded from candidate scoring and selection, and the selected formulation was subsequently evaluated against the held-out study. (a) Scores of V1–V6 under the primary endpoint-to-curve weighting of 75:25; boxes denote the selected lowest-scoring formulations. (b) Formulations selected using endpoint-to-curve error weights of 75:25, 50:50, and 25:75. (c) Held-out endpoint MAPE and curve normalized root-mean-square error (NRMSE) of the selected formulation. Under the primary weighting, V6 was reselected in four of the five folds (Figure 3(a)). The same four folds selected V6 under all three weighting schemes (Figure 3(b)), and their held-out endpoint MAPEs ranged from 3.2% to 13.9% (Figure 3(c)). Selection changed only when the Tjelta study [41] was excluded. In this fold, the 75:25 and 50:50 weightings selected V4, whereas the 25:75 weighting, which placed greater emphasis on the process curves, selected V5 (Figure 3(b)). Tjelta provides the only full-scale field heave–depth profile in the evidence pool. Its exclusion left the candidate scores primarily determined by centrifuge and laboratory model tests, which favored V4 under endpoint-dominated evaluation and V5 under curve-dominated evaluation. When evaluated against the held-out Tjelta profile, the selected V4 produced an endpoint error of 45.5% and a curve normalized root-mean-square error (NRMSE) of 39.2%. For comparison, V6 produced corresponding errors of 8.2% and 20.4%, respectively. Overall, candidate selection remained stable across four of the five study exclusions and all tested weighting schemes within those four folds, but changed when the only full-scale study was removed. The relatively lower held-out errors of V6 for the Tjelta profile indicate that it provided a more consistent balance across centrifuge tests, laboratory model tests, and the available full-scale field response, supporting its retention as the selected formulation. The selected model integrates three plug-heave mechanisms The selected V6 formulation predicts plug heave by combining skirt-penetration-induced geometric displacement (MG 29005_ 28999), seepage-driven soil void-ratio evolution (MS 29005_ 29011), and near-wall soil dilation (MD 29005_ 28996) (Figure 4). Figure 4: Physical structure and coordinate definitions of the selected soil-plug formulation. The fixed coordinate y 29049 is measured downward from the initial seabed, and the local coordinate x 29048 is measured downward from the current plug surface. At penetration depth z 29050, the plug length is Lplug=z+hea|e 29004_ 29040 29036 29045 29031 12349 29050 8235 29032_ 29032 29029 29025 69640972 29029. The initial material coordinate ζ0 28944_ 28720 tracks the soil entering the internal control volume. MG 29005_ 28999, MS 29005_ 29011, and MD 29005_ 28996 denote skirt-penetration-induced geometric displacement, seepage-driven soil void-ratio evolution, and near-wall soil dilation, respectively. Dashed blue curves schematically indicate the seepage paths. At penetration depth z 29050, these three contributions are assembled through the solid-volume-conservation relation: hea|e(z)=∫0z[αA1+e(ζ0,z)1+e0−1]dζ0, 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 12349 4946 _ 28720 29050 67482370 28939_ 28993 28721 8235 29029 67273472 28944_ 28720 24891 29050 84054785 28721 8235 29029_ 28720 8704 28721 84267779\, 29028 28944_ 28720 24891 (1) where ζ0 28944_ 28720 labels a soil layer by its initial depth below the seabed, and z 29050 denotes the current penetration depth. The quantity e(ζ0,z) 29029 67273472 28944_ 28720 24891 29050 84054785 is the void ratio of that soil layer after the MS 29005_ 29011 update and subsequent MD 29005_ 28996 correction at penetration depth z 29050. The initial void ratio is e0 29029_ 28720, and αA=(Do/Di)2 28939_ 28993 12349 67273472 28996_ 29039 68408078 28996_ 29033 84054785 28722 is the caisson area ratio, where Do 28996_ 29039 and Di 28996_ 29033 are the outer and inner caisson diameters, respectively. When no void-ratio change occurs, e(ζ0,z)=e0 29029 67273472 28944_ 28720 24891 29050 84054785 12349 29029_ 28720, and Equation (1) reduces to the geometric limit: hea|e(z)=(αA−1)z, 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 12349 67273472 28939_ 28993 8704 28721 84054785 29050 24891 (2) Equation (2) recovers the geometric contribution MG 29005_ 28999. In the full formulation, MS 29005_ 29011 and MD 29005_ 28996 modify the predicted plug heave through their updates to e(ζ0,z) 29029 67273472 28944_ 28720 24891 29050 84054785, while Equation (1) assembles these contributions into the observable response in V6. The coupled seepage, effective-stress, and void-ratio relations used to obtain e(ζ0,z) 29029 67273472 28944_ 28720 24891 29050 84054785 are provided in Supplementary Note 1. Supplementary Fig. 1 presents the sequence through which these component calculations are evaluated and assembled at each penetration step. Multi-agent audits recover known and previously unflagged issues The blinded replay tested whether the retained development records contained enough information for an isolated audit role to recover known modeling issues. With the original discussions and expected answers withheld, the audit role recovered and localized all 9 issues, reported no issues outside the target set, and completed the scored replay without human assistance (Table 3). These results show that the retained records were sufficient to identify where each known issue arose and which physical rule, computational operation, or code segment required re-examination. The derivation-to-code probes then examined whether a separate audit role could detect implementation problems missed by the predefined checks. Across the 5 probes, the implementation role passed all 36 fixed checks. The audit role then reported five additional findings across 4 of the 5 probes (Table 3). These findings concerned an omitted rule for refreshing the final output, branch execution, the scope of a stated rule, an incomplete audit trace, and a numerical safeguard applied to the wrong term in an equation. This indicates that the role-separated audit extended the evaluation beyond the coverage of the predefined checks, recovering an omitted requirement and identifying implementation issues that those checks did not address. After these findings had been addressed, the final code used for benchmark evaluation passed all 37 code tests, comprising 33 module-level tests and 4 regression tests (Table 3). It also passed all 5 independently calculated reference checks. Within the scope of these tests, the results confirmed that the final implementation followed the equations and calculation sequence described in the manuscript. The blinded-replay and derivation-to-code probe records are provided in Supplementary Note 2, and the final-code tests and regression baselines are reported in Supplementary Note 3. Table 3: Results of workflow audit and final code verification. Audit task Observed result Supported conclusion Blinded replay of known issues All 9 known issues were recovered and localized; no issues outside the target set, and no human assistance was provided during the scored replay. The retained development records enabled an isolated audit role to recover and localize known issues under blinded conditions. Derivation-to-code probe audit All 36 predefined checks passed; the separate audit reported 5 additional findings across 4 of the 5 probes, with no finding in one probe. The role-separated audit identified both an omitted predefined rule and implementation issues not covered by the predefined checks. Final code verification All 37 code tests passed, comprising 33 module-level tests and 4 regression tests; all 5 independently calculated reference checks also passed. Within the checked scope, the final implementation followed the equations and calculation sequence described in the manuscript. Staged benchmarks support final-state and process-level predictions We evaluated V6 using three levels of plug-heave evaluation: formula reduction and volume closure, final state prediction, and continuous process prediction (Figure 5(a)). Figure 5: Staged benchmark evaluation of V6 plug-heave predictions. (a) Evidence map linking the available source data, benchmark problem, and supported claim at each evaluation level. (b) Formula reduction and volume closure for 7 cases using source-reported changes in void ratio or volumetric strain. The diagonal line represents equality between predicted and measured final heave, and the shaded region denotes the predefined ±10% comparison band. (c) Final state prediction for 9 benchmark cases. Black diamonds are measured final heave, blue points represent central predictions, and horizontal bands show the sensitivity range obtained by varying eUNKNOWN/emax 29029_ 28672 68408078 29029_ 29037 29025 29048 from 0.9 to 1.1. Percentages report the signed errors of the central predictions. (d) Predicted and measured plug-heave–depth histories for 5 cases, together with the corresponding sensitivity bands and final-point errors. The bands in panels (c) and (d) represent prescribed parameter sensitivity rather than probabilistic confidence intervals. At the first level, final plug heave was calculated using reported changes in void ratio or volumetric strain. 6 of the 7 V6 calculations fell within or close to the predefined ±10% 8710 28721 28720\% comparison band (Figure 5(b)). 6 cases used void-ratio changes reported by Koteras and Ibsen [42], while the Ragni et al. [15] T3 case used an image-derived volumetric strain. Only T4 showed a clear deviation, with an error of +29.5% 8235 28722 28729 314 28725\%. This difference may reflect experimental scatter or uncertainty in the reported measurement: T4 and T5 had the same void-ratio change and similar penetration depths, but their measured final plug heaves differed by 7 m, equivalent to 16.7% of the T4 measurement. The Ragni T3 case yielded an error of +6.6% 8235 28726 314 28726\%, showing that the volume-conservation relation also reproduced the measured final heave when volumetric strain was provided directly. At the second level, V6 predicted final plug heave for each case from the initial soil state and suction-pressure history. The prediction sensitivity bands encompassed 7 of the 9 benchmark endpoints, and the central predictions yielded a MAPE of 12.4% 28721 28722 314 28724\% (Figure 5(c)). The endpoints for Kim and Kim [16] T1-1 S 25g and Tran and Randolph [14] TRAN-T3 lay outside these bands. The upper prediction bounds from the predefined sensitivity analysis reduced the respective shortfalls from 26.9% 28722 28726 314 28729\% to 11.7% 28721 28721 314 28727\% and from 28.1% 28722 28728 314 28721\% to 17.1% 28721 28727 314 28721\%, respectively. The reference endpoint for TRAN-T3 was inferred from the final digitized penetration depth under a plug-filled endpoint assumption, which introduced additional uncertainty into the endpoint comparison. Overall, V6 reproduced final plug heave for most benchmark cases, while the sensitivity analysis reduced the underprediction in the two outside-band cases. The third level extended the prediction from penetration endpoints to complete plug-heave–depth histories. Across 5 cases, V6 captured the overall increase in plug heave with penetration, and all 5 measured endpoints lay within the predefined prediction sensitivity bands (Figure 5(d)). Final-point errors ranged from 1.3% 28721 314 28723\% to 8.2% 28728 314 28722\%, giving a MAPE of 4.2% 28724 314 28722\%. V6 agreed most closely with the three curves reported by Zhao et al. [43], with curve normalized root-mean-square errors (NRMSEs) of 3.9% 28723 314 28729\%–4.3% 28724 314 28723\%. The Wang et al. [44] and Tjelta [41] cases yielded curve NRMSEs of 10.8% 28721 28720 314 28728\% and 20.4% 28722 28720 314 28724\%, respectively. The largest process-level discrepancy occurred during early penetration in the Tjelta case, for which V6 predicted a faster increase in plug heave than was measured. This early overprediction is consistent with the treatment of geometric displacement in V6: soil displaced by skirt penetration contributes immediately to the calculated plug heave, while under the relatively low suction present during early penetration, some displaced soil may first redistribute before producing measurable heave. As penetration continued, the predicted and measured responses converged and reached comparable final magnitudes. Together, the three benchmark levels support distinct aspects of V6 rather than a single undifferentiated validation claim. The first evaluates volume closure when the end-state soil response is supplied, the second evaluates final-heave prediction from the initial state and suction history, and the third evaluates the evolution and final magnitude of plug heave. Detailed case results are reported in Supplementary Table 13 for formula reduction and volume closure, Supplementary Table 14 for final state prediction, and Supplementary Tables 15 and 16 for continuous process prediction. Discussion We used a multi-agent workflow operating within human-defined physical boundaries to connect source-supported knowledge, theory construction, numerical implementation, audit, and benchmark evaluation. Iteration through 6 candidate formulations produced V6, which integrates skirt-penetration-induced geometric displacement, seepage-driven soil void-ratio evolution, and near-wall soil dilation through solid-volume conservation. Its central predictions yielded a MAPE of 12.4% across 9 final-state benchmark cases and a final-point MAPE of 4.2% across 5 measured plug-heave–penetration-depth histories. With these, this study produced two linked outcomes: a physically interpretable model of soil-plug evolution and a traceable route from its knowledge sources to its numerical implementation and evaluation. This study extends agent-assisted scientific computing from isolated retrieval, coding, or data-analysis tasks to a longer engineering-modeling chain. Existing multi-agent systems have generated computational strategies within predefined scientific machine-learning problems or linked literature-grounded hypotheses to experimental data analysis [45, 46]. In the present study, no complete soil-plug model was available at the outset. In this case, the workflow had to identify relevant physical mechanisms, derive their governing relations, discretize the equations, and implement the resulting solver. Human experts retained authority over admissible knowledge sources, modeling assumptions, and model revisions. When benchmark evaluation revealed a substantial discrepancy between predicted and measured plug heave, or a trade-off between final-heave and process errors, the agents formulated the next physical question and retrieved relevant source-supported relations. Each admitted relation was then translated into equations, numerical operations, and code implementation. The methodological contribution therefore lies not in removing human judgment, but in organizing human and agent roles so that the complete modeling process remains inspectable. The evolution from V1 to V6 illustrates how predictive discrepancies can be converted into physically meaningful diagnostic questions. V1–V3 continued to underpredict plug heave after seepage and confined unloading had been introduced, indicating that their stress updates generated insufficient void-ratio change. The workflow consequently examined how movement along the stress path should be converted into volumetric response. Introducing a state-dependent volume-change relation in V4 substantially reduced both final-heave and process final-point errors. V5 further reduced the process final-point error but increased the final-heave error, revealing a trade-off between the two evaluation targets. V6 addressed this trade-off by integrating the admitted seepage, soil-state, and near-wall dilation mechanisms while maintaining consistent final-state updating. In this way, each model revision was linked to a defined physical question, supporting evidence, a mathematical expression, and an implementation path. The resulting formulation also provides a physical interpretation of plug-heave development. Skirt penetration determines the volume of soil entering the caisson, while the suction-induced stress path changes the soil state and near-wall dilation further modifies the void-ratio field. Solid-volume conservation then converts this evolving field into observable plug heave. Therefore, the plug heave can be interpreted as a process response initiated by geometric displacement and modulated by suction-driven soil-state evolution and interface behavior. Complementary benchmark evidence supported the selection of V6. Its prediction sensitivity bands encompassed 7 of the 9 final-state benchmark cases and the endpoints of all 5 measured plug-heave–penetration-depth histories. The retrospective leave-one-study-out analysis reselected V6 in 4 of the 5 folds and under all 3 scoring schemes in those folds. Selection changed only when the full-scale Tjelta study [41] was excluded. The formulation selected without this field record produced substantially larger errors than V6 when evaluated against the held-out Tjelta profile. This comparison highlights the importance of full-scale process data in constraining model selection and supports the balance achieved by V6 across final-state, process-level, and experimental-scale evidence. From an engineering perspective, V6 extends the conventional geometric heave estimate to a process-level prediction controlled jointly by the suction-pressure history and evolving soil state. It retains a transparent geometric baseline while describing how plug heave develops throughout penetration. The resulting formulation provides a physically interpretable reduced-order model for assessing both final plug heave and its development during suction-caisson installation. Independent role-separated auditing addresses a different question from benchmark evaluation: whether the declared physical relations were translated consistently into equations and executable code. The blinded replay recovered and localized all 9 known problems without access to the expected answers, showing that the retained records preserved the information required for subsequent examination. The derivation-to-code audit then identified 5 additional findings across 4 probes after all 36 predefined implementation checks had passed. This indicates that the separate audit role extended the coverage of the predefined checks by identifying an omitted requirement and cross-layer inconsistencies in the transition from derivation to implementation. Benchmark evaluation examined what the solver predicted, while the audit examined how the solver was constructed. The evidence boundary of the current study is determined primarily by the available process data and the scope of the demonstration. Only 5 comparable plug-heave–penetration-depth histories are available, including one full-scale field record, and few studies jointly report suction-pressure history, continuous plug heave, and independently measured soil state. V6 currently applies to vertical penetration in saturated homogeneous sand and assumes pure soil inflow, quasi-steady Darcy seepage, and cross-sectionally averaged near-wall soil dilation. More complete field observations and explicit treatment of early soil redistribution, stratification, and three-dimensional interface behavior will further test and extend both the physical model and the workflow. In summary, this study treats engineering-model development as an inspectable scientific process rather than an unreported precursor to the final solver. The selected formulation (V6) demonstrates that this workflow can produce an engineering model with explicit physical meaning, while benchmark evaluation and independent audit examine its predictive outcome and formation process, respectively. The central implication is that AI agents can participate in constructing engineering theory within explicit human-defined physical boundaries, while the model’s knowledge sources, derivation path, numerical implementation, and predictive performance remain open to scientific scrutiny. Methods Human-guided physics-constrained multi-agent workflow The human-guided physics-constrained multi-agent workflow organizes model development into five linked stages: knowledge admission, skill composition, solver construction, independent audit, and benchmark evaluation (Fig. 6). These stages are implemented through the knowledge engine, constraint and skill-composition modules, execution layer, audit layer, and benchmark-validation interface. Figure 6: Overview of the human-guided physics-constrained multi-agent workflow. Blue and orange dashed arrows indicate feedback on physical assumptions and implementation issues, respectively. The workflow is designed to identify and avoid errors that can enter the modeling chain before the solver is validated against benchmark problems. These errors mainly arise from three aspects: 1. What knowledge is admitted, e.g., an assumption used outside its admissible range. 2. How agents use the knowledge, e.g., a symbol that changes meaning during derivation, or a temporary iteration value recorded as a final input. 3. How solver outputs are verified, e.g., a failure to trace a reported prediction back to the physical relations and calculation steps that produced it. Figure 6 maps these critical aspects onto the workflow structure: the knowledge engine supplies traceable physical records; the constraint layer decides which physical laws and operations are admissible; the execution layer turns the admitted components into equations and code under fixed notation and state rules; and the audit layer checks consistency across the physical rules, equations, and code and returns unresolved inconsistencies for human review. Throughout this sequence, human experts retain authority over record admission and model revision, while agents perform bounded retrieval, derivation, code implementation, and auditing tasks. The first human decision occurs in the knowledge engine, the entry point for physical information. It contains two categories of agentic skills that serve as building blocks for constructing the theoretical framework. 1. Static skills: human-approved physical laws, constitutive relationships, or numerical operators. An agent can add a knowledge block to the model only when its assumptions, variables, and boundary conditions are compatible with those of its neighboring blocks. This bounded search-and-use strategy forms the skill-RAG operation of the workflow. The language model is guided by records retrieved during use [23], and the retrieval target is a task-relevant skill rather than unrestricted text [24]. 2. Dynamic retrieval: if the initially defined static skills do not contain everything needed for model assembly, an agent may use the LLM’s pretrained knowledge, web search, or other tools to locate candidate theories, published mechanisms, closure assumptions, or numerical operators [20]. Before use, the LLM-identified item needs to be reframed as a compatible skill block. This process requires its source to be identified, its physical relation to be defined, and its allowable use to be stated, including variables, units, assumptions, and applicable range. A human reviewer then decides whether the item is admitted as a new skill or kept only as supporting evidence. After records have been admitted, the execution layer controls how agents use them to build the solver. Agents select compatible records, assemble mechanisms, derive equations, and implement code modules under the approved assumptions and parameter rules. The execution layer also fixes numerical conventions before solver outputs are accepted, including coordinates, signs, gradient directions, and the volumetric-strain convention. It separates trial values from the accepted state, defined here as the checked version of variables, inputs, and update history submitted for reporting. This rule ensures that version records, committed histories, and reported outputs refer to the same calculation iteration. The audit layer then evaluates whether the continuous equations, discrete update rules, and code modules remain consistent. If a mismatch is found, the affected item returns to the execution layer for correction and is then checked again. Agents can identify and localize these inconsistencies, but they cannot admit new records or authorize model changes until human review. In this implementation, the derivation, coding, and audit roles were run in OpenAI Codex with GPT-5.4 as the underlying model. Codex performed bounded tasks, while humans decided whether new records could be admitted and whether the model should be revised. The retained audit materials and results are documented in Supplementary Note 2. The same loop applies after benchmark validation. After the audit layer establishes internal consistency, the model parameters remain fixed while differences between predicted and reported plug heave are evaluated against two predeclared bounds: the comparison band, which defines acceptable differences from the literature data, and the prediction sensitivity band, which represents prediction variation caused by documented input uncertainty. An outside-bound prediction attributed to a missing mechanism or unsuitable assumption returns the model to the knowledge engine, constraint layer, and human review. A difference attributable to documented source uncertainty, benchmark definition, or declared model scope is retained as a bounded modeling limitation. Soil-plug model-development task To apply the workflow, we defined soil-plug evolution during vertical suction-caisson penetration in saturated, homogeneous sand as the model-development task (Fig. 1(b)). The modeled system was the soil entering the caisson and forming an internal plug during penetration, and the target response was its rise relative to the surrounding seabed. Task inputs comprised caisson geometry, initial soil state and material properties, penetration depth, and the available suction-pressure history. From these inputs, the agents were tasked with constructing an executable model to calculate plug heave as penetration progressed and at the final penetration depth. At task initiation, the admitted static skills comprised solid-volume conservation, Darcy-type seepage, critical-state-line (CSL) descriptions, and Bolton-type dilatancy. These skills formed the initial knowledge base for model development. When additional physical relations or closure assumptions were required during model development, they followed the knowledge-admission procedure described in the preceding subsection. Development time was recorded for the derivation, implementation, and independent audit roles after the physical knowledge and model inputs had been prepared. For each version, role duration was the elapsed time required to complete the corresponding output, and total agent execution time was the sum of the three role durations. Workflow audit and implementation verification To evaluate whether the workflow could detect and localize problems along the model-development chain, we used two audit tasks followed by a final code check: 1. Blinded audit replay: tested whether an audit role could recover known issues using only the retained development materials. 2. Derivation-to-code probes: tested whether a separate audit role could identify problems not covered by the fixed checks after code generation. 3. Final code verification: checked whether the corrected code used for benchmark evaluation followed the equations and calculation sequence described in the manuscript. In the blinded audit replay, the audit role received development materials but was not given the target issue type, location, or answer key. The task tested whether the role could recover and localize the issue from those materials alone. We first converted 9 consistency issues identified during workflow execution into audit packets that contained no answers, while retaining the corresponding answer keys separately for scoring. Each packet contained only the relevant physical relation, computational rule, and pre-correction code. An isolated audit role then reviewed one packet using a general audit instruction, without access to issue labels, answer keys, prior discussions, or any material outside the packet. Finally, we compared each audit report with the corresponding answer key and scored whether the role detected and localized the issue, explained its consequences, and proposed a valid correction or follow-up check. False positives and human assistance were recorded separately. The derivation-to-code probes extended the assessment to code generation. We constructed 5 controlled probes. For each probe, an implementation role generated code from the supplied equations and state-update rules, followed by the corresponding fixed checks. Across the 5 probes, 36 predeclared checks were completed in total. A separate audit role then received the derivation materials, generated code, and check records. It assessed equation-to-code consistency, focusing on boundary and branch handling, the timing of state and output updates, and the completeness of the check records. Full procedures for preparing the audit materials, limiting the information provided to each role, and scoring the reports are described in Sections S2.1–S2.2 of Supplementary Note 2. After the audit findings had been addressed, the code used for benchmark evaluation was checked for consistency with the governing equations and calculation sequence. The check covered key calculations, repeated runs with fixed inputs, and comparisons between key outputs and independently calculated reference values. The complete test list and passing criteria are provided in Supplementary Note 3. Benchmark evaluation and model selection Benchmark datasets For the soil-plug case study, 3 benchmark problems were formulated according to the type of suction-response information available in the literature: reported end-state volumes, final plug heave, or complete heave–depth histories. 1. Formula reduction and volume closure were evaluated using 7 cases from 2 studies [42, 15]. Void ratios or volumetric strains before and after suction installation were provided to calculate the final plug heave directly from the end-state volumes and to verify the overall volume closure. 2. Final state prediction was evaluated using 9 cases from 2 studies [16, 14]. The pre-suction soil void ratio, the suction-pressure history, and the test conditions were provided to predict the final plug heave and the associated prediction sensitivity band. 3. Continuous-process prediction was evaluated using 5 cases with measured heave–depth histories from 3 studies [41, 44, 43]. The initial soil state and the suction-pressure history were provided to predict the development of plug heave during penetration. Input and parameter assignment Benchmark inputs were assigned through a fixed source hierarchy. Validation targets were obtained from reported measurements or reconstructed from published figures. Values explicitly reported in the source studies were used first: measured quantities (direct) and test-design parameters (design parameter). Inputs not reported explicitly were derived where possible from reported information or physical relationships (derived). The remaining inputs were assigned from data for the same material (analog) or from predefined engineering assumptions (assumption). Each adopted value and its source class were recorded in a parameter-source ledger (Supplementary Data 2). The label unavailable identifies entries for which no sourced value was adopted; the label mixed identifies rows combining values from multiple source classes. All inputs were fixed before any model calculation, so that no input could be adjusted in response to prediction outcomes. Each candidate formulation was then evaluated by comparing its predictions with the validation targets. Supplementary Note 4 documents the inputs, source labels, and derivations for each benchmark case, and the underlying numerical records are provided as Supplementary Data 1. Most benchmark inputs were taken directly from the source studies or derived through the hierarchy above. Three model parameters, however, were assigned by common rules in the final-state and continuous-process predictions. None of the five prediction studies provided case-specific values of the critical-state-line intercept ratio eUNKNOWN/emax 29029_ 28672 68408078 29029_ or Poisson’s ratio ν 28951; the central value of eUNKNOWN/emax 29029_ 28672 68408078 29029_ was therefore set to 1.0, the midpoint of the sensitivity range adopted below, and ν 28951 was set to 0.30 [37, 47]. Kim and Kim [16] reported case-specific ki/kout 29035_ 29033 68408078 29035_ 29039 29045 29044 values ranging from approximately 1.26 to 1.46, and a representative value of 1.3 was adopted for their cases; the remaining cases used ki/kout=3.0 29035_ 29033 68408078 29035_ 29039 29045 29044 12349 28723 314 28720, following Houlsby and Byrne [6]. Formula reduction and volume closure required no such assignments, because their calculations used the source-reported end-state volumes directly. Sensitivity to eUNKNOWN/emax 29029_ 28672 68408078 29029_ was evaluated by varying it by ±10% 8710 28721 28720\% about the central value, giving a range of 0.9–1.1; the adopted range was based on reported critical-state-line and state-parameter data for sand [48, 49, 50, 51]. The resulting spread defines the prediction sensitivity band for the final-state and continuous-process predictions. This band isolates the effect of eUNKNOWN/emax 29029_ 28672 68408078 29029_ : no case-specific values of this ratio were available for the prediction studies, and it sets the low-stress reference used in the void-ratio update, whereas Poisson’s ratio and the permeability ratio were held at the values assigned above. Evaluation metrics The metrics were chosen to match the validation data available for each benchmark. For formula reduction and volume closure, final plug heave was calculated from the reported void-ratio or volumetric-strain information and compared with the reported final plug heave. Agreement was assessed against the predefined comparison band of ±10% 8710 28721 28720\%. For final state prediction, we expressed the signed percentage error for benchmark case c 29027 as Ec=hpred,c−hmeas,chmeas,c×100%, 28997_ 29027 12349 29032_ 29040 29042 29029 29028 24891 29027 8704 29032_ 29037 29029 29025 29043 24891 29027 29032_ 29037 29029 29025 29043 24891 29027 8706 28721 28720 28720\% 24891 (3) where hpred,c 29032_ 29040 29042 29029 29028 24891 29027 is the predicted final plug heave and hmeas,c 29032_ 29037 29029 29025 29043 24891 29027 is the corresponding final plug-heave value used for validation. The value of hmeas,c 29032_ 29037 29029 29025 29043 24891 29027 was either reported directly or reconstructed from the final digitized penetration depth for that case, as described above. A positive Ec 28997_ 29027 indicates overprediction, whereas a negative value indicates underprediction. We then summarized the absolute case-level errors as the mean absolute percentage error: MAPE=1Nc∑c=1Nc|Ec|, 29005 28993 29008 28997 12349 28721 29006_ 29027 4944 _ 29027 12349 28721 29006_ 29027 69640972 28997_ 29027 69640972 24891 (4) where Nc 29006_ 29027 is the number of cases included in the final state prediction average. For continuous process prediction, the predicted plug heave was evaluated at the measured penetration depths. The error at the final measured depth was computed using Equation (3). Averaging its absolute value across the process cases gives the final-point MAPE. Because percentage error is undefined when the measured heave is zero, curve MAPE was calculated only at points with nonzero measured plug heave: MAPEcur|e=1N+∑p∈+|hpred,p−hmeas,phmeas,p|×100%, 29005 28993 29008 28997_ 29027 29045 29042 69640972 29029 12349 28721 29006_ 8235 4944 _ 29040 12850 29008_ 8235 69640972 29032_ 29040 29042 29029 29028 24891 29040 8704 29032_ 29037 29029 29025 29043 24891 29040 29032_ 29037 29029 29025 29043 24891 29040 69640972 8706 28721 28720 28720\% 24891 (5) where p 29040 indexes the process points, and hpred,p 29032_ 29040 29042 29029 29028 24891 29040 and hmeas,p 29032_ 29037 29029 29025 29043 24891 29040 are the predicted and measured plug heave at point p 29040, respectively. The set + 29008_ 8235 contains the points for which both values are finite and hmeas,p/=0 29032_ 29037 29029 29025 29043 24891 29040 566 61 28720, and N+=|+| 29006_ 8235 12349 69640972 29008_ 8235 69640972. To complement this relative measure, we calculated the root-mean-square error in physical units over all process points with finite predicted and measured plug heave, including points with zero measured heave: RMSEcur|e=[1Np∑p∈f(hpred,p−hmeas,p)2]1/2, 29010 29005 29011 28997_ 29027 29045 29042 69640972 29029 12349 67482370 28721 29006_ 29040 4944 _ 29040 12850 29008_ 29030 67273472 29032_ 29040 29042 29029 29028 24891 29040 8704 29032_ 29037 29029 29025 29043 24891 29040 84054785 28722 84267779 28721 68408078 28722 24891 (6) where f 29008_ 29030 is the set of process points for which both the predicted and measured plug heave are finite, and Np=|f| 29006_ 29040 12349 69640972 29008_ 29030 69640972. Plug heave is expressed in centimeters. Candidate screening Each candidate formulation generated by the workflow was evaluated using the same benchmark cases, input assignments, and parameter rules. Because the candidate formulations shared the same volume-closure relation, comparison during model development focused on final state prediction and continuous process prediction. Final-heave MAPE characterized agreement with the endpoint data, whereas process final-point MAPE characterized agreement at the end of the measured penetration histories. Both metrics were retained in the development records and used to determine whether the formulation entered another cycle of audit and revision. LOSO selection stability The dataset used for candidate selection contained 5 source studies. Each study was held out once, giving 5 retrospective leave-one-study-out (LOSO) folds. In each fold, the remaining 4 studies were used to rescore the candidate models and identify the best-performing model, which was then evaluated against the held-out study. A change in the selected model after excluding a study indicated sensitivity to that study, whereas consistent selection across the 5 folds indicated stable model selection. Candidate scores combined final and process performance. The final-heave component, Efinal 28997_ 29030 29033 29038 29025 29036, was calculated from the absolute percentage errors in final plug heave. The process component, Ecur|e 28997_ 29027 29045 29042 69640972 29029, measured the error over the complete heave–depth curve and included only studies that reported process histories. To prevent studies with more cases or curves from dominating the score, errors were first averaged within each study and then across studies. For LOSO scoring, each curve RMSE was normalized by the measured final plug heave: NRMSEcur|e=100RMSEcur|ehmeas,final. 29006 29010 29005 29011 28997_ 29027 29045 29042 69640972 29029 12349 28721 28720 28720 29010 29005 29011 28997_ 29027 29045 29042 69640972 29029 29032_ 29037 29029 29025 29043 24891 29030 29033 29038 29025 29036 314 (7) where RMSEcur|e 29010 29005 29011 28997_ 29027 29045 29042 69640972 29029 is the root-mean-square difference between the predicted and measured curves, and hmeas,final 29032_ 29037 29029 29025 29043 24891 29030 29033 29038 29025 29036 is the measured plug heave at the end of the curve. Averaging these normalized errors within each study and then across the contributing studies gave Ecur|e 28997_ 29027 29045 29042 69640972 29029. The primary score assigned 75% of the weight to final-heave error and 25% to process-curve error: J=0.75Efinal+0.25Ecur|e. 29002 12349 28720 314 28727 28725 28997_ 29030 29033 29038 29025 29036 8235 28720 314 28722 28725 28997_ 29027 29045 29042 69640972 29029 314 (8) The candidate model with the lowest J 29002 was selected in each fold. The analysis was repeated using final-heave and process-curve weights of 50:50 and 25:75 to assess whether model selection depended on the relative weighting of the two types of evidence. Data availability The published experimental data and model predictions used in the benchmark validation are provided as Supplementary Data 1. The model inputs and their sources are provided as Supplementary Data 2. Code availability The code used to run the final soil-plug solver and reproduce the benchmark predictions is publicly available under the MIT License at https://github.com/shijie774302906/soil-plug-solver. References [1] International Energy Agency, Offshore wind outlook 2019, Tech. rep., International Energy Agency, Paris (2019). URL https://w.iea.org/reports/offshore-wind-outlook-2019 [2] N. Gonzalez, P. Serna-Torre, P. A. S’anchez-P’erez, R. Davidson, B. Murray, M. Staadecker, J. Szinai, R. Wei, D. M. Kammen, D. A. Sunter, P. Hidalgo-Gonzalez, Offshore wind and wave energy can reduce total installed capacity required in zero-emissions grids, Nature Communications 15 (1) (2024). doi:10.1038/s41467-024-50040-6. [3] Y. Zhao, Y. Tao, Y. Chen, J. Yan, Z. Zeng, et al., Increasing extreme winds challenge offshore wind energy resilience, Nature Communications 16 (2025) 9529. doi:10.1038/s41467-025-65105-3. [4] X. Wu, et al., Foundations of offshore wind turbines: A review, Renewable and Sustainable Energy Reviews 104 (2019) 379–393. doi:10.1016/j.rser.2019.01.012. [5] OWA, Suction installed caisson foundations for offshore wind: Owa design guidelines, Tech. Rep. Issue 1.0, Offshore Wind Accelerator, The Carbon Trust (2019). [6] G. T. Houlsby, B. W. Byrne, Design procedures for installation of suction caissons in sand, Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 158 (3) (2005) 135–144. doi:10.1680/geng.158.3.135.66297. [7] O. Harireche, M. T. Naqash, Q. U. Farooq, A full numerical model for the installation analysis of suction caissons in sand, Ocean Engineering 234 (2021) 109173. doi:10.1016/j.oceaneng.2021.109173. [8] O. C. Zienkiewicz, A. H. C. Chan, M. Pastor, B. A. Schrefler, T. Shiomi, Computational Geomechanics with Special Reference to Earthquake Engineering, John Wiley & Sons, Chichester, 1999. [9] A. Gens, Soil–environment interactions in geotechnical engineering, Géotechnique 60 (1) (2010) 3–74. doi:10.1680/geot.9.P.109. [10] G. Wilson, D. A. Aruliah, C. T. Brown, N. P. Chue Hong, M. Davis, R. T. Guy, S. H. D. Haddock, K. D. Huff, I. M. Mitchell, M. D. Plumbley, B. Waugh, E. P. White, P. Wilson, Best practices for scientific computing, PLOS Biology 12 (1) (2014) e1001745. doi:10.1371/journal.pbio.1001745. [11] H. C. Ates, P. Q. Nguyen, L. Gonzalez-Macia, E. Morales-Narv’aez, F. G"uder, J. J. Collins, C. Dincer, End-to-end design of wearable sensors, Nature Reviews Materials 7 (2022) 887–907. doi:10.1038/s41578-022-00460-x. [12] Z. Hu, W. Yan, Data-driven modeling of process-structure-property relationships in metal additive manufacturing, npj Advanced Manufacturing 1 (2024) 3. doi:10.1038/s44334-024-00003-y. [13] I. T. Song, J. Kang, J. Koh, H. Choi, H. Yang, E. Park, J. Lee, W. Cho, Y.-m. Lee, S. Lee, et al., Thermal runaway prevention through scalable fabrication of safety reinforced layer in practical li-ion batteries, Nature Communications 15 (2024) 8294. doi:10.1038/s41467-024-52766-9. [14] M. N. Tran, M. F. Randolph, Variation of suction pressure during caisson installation in sand, Géotechnique 58 (1) (2008) 1–11. doi:10.1680/geot.2008.58.1.1. [15] R. Ragni, B. Bienen, S. Stanier, C. O’Loughlin, M. Cassidy, Observations during suction bucket installation in sand, International Journal of Physical Modelling in Geotechnics 20 (3) (2020) 132–149. doi:10.1680/jphmg.18.00071. [16] J. H. Kim, D. S. Kim, Soil plug heave induced by suction bucket installation on sand via centrifuge model tests, Marine Georesources & Geotechnology 38 (10) (2020) 1245–1256. doi:10.1080/1064119X.2019.1675204. [17] J. Shi, B. A. Malik, Z. Liang, L. Jing, Y. Lu, F. Gao, Z. Ouyang, Installation effects on the performance of hybrid suction caissons in sand, Ocean Engineering 359 (2026) 126025. [18] P. Zhang, L. Mu, Y. Lu, M. Huang, J. Sun, Microscopic insights into suction bucket installation in sand: Coupled coarse-grained cfd-dem simulations, Computers and Geotechnics 167 (2024) 106060. [19] L. Mu, T. Zhou, Y. Lu, G. Yu, J. Sun, Experimental investigation of suction-assisted penetration of compartmented bucket foundations, Applied Ocean Research 166 (2026) 104923. [20] S. Yao, J. Zhao, D. Yu, N. Du, I. Shafran, K. Narasimhan, Y. Cao, React: Synergizing reasoning and acting in language models, proceedings of ICLR 2023. arXiv:2210.03629 (2023). [21] D. Guo, D. Yang, H. Zhang, et al., Deepseek-r1 incentivizes reasoning in llms through reinforcement learning, Nature 645 (8081) (2025) 633–638. doi:10.1038/s41586-025-09422-z. [22] D. A. Boiko, R. MacKnight, B. Kline, G. Gomes, Autonomous chemical research with large language models, Nature 624 (7992) (2023) 570–578. doi:10.1038/s41586-023-06792-0. [23] P. Lewis, E. Perez, A. Piktus, F. Petroni, V. Karpukhin, N. Goyal, H. Kuttler, M. Lewis, W.-t. Yih, T. Rocktaschel, S. Riedel, D. Kiela, Retrieval-augmented generation for knowledge-intensive nlp tasks, in: Advances in Neural Information Processing Systems, Vol. 33, 2020, p. 9459–9474, arXiv:2005.11401. [24] K. Wei, R. Li, X. Zhu, Z. Xue, J. Han, J. Niu, F. Yang, Skill-rag: Failure-state-aware retrieval augmentation via hidden-state probing and skill routing, arXiv:2604.15771 (2026). doi:10.48550/arXiv.2604.15771. [25] S. Li, T. Marwah, J. Shen, W. Sun, A. Risteski, Y. Yang, A. Talwalkar, Codepde: An inference framework for llm-driven pde solver generation, Transactions on Machine Learning Research. arXiv:2505.08783 (2026). [26] T. Kim, T. S. Yun, H. S. Suh, Can chatgpt implement finite element models for geotechnical engineering applications?, International Journal for Numerical and Analytical Methods in Geomechanics 49 (6) (2025) 1747–1766. doi:10.1002/nag.3956. [27] A. M. Bran, S. Cox, O. Schilter, C. Baldassari, A. D. White, P. Schwaller, Augmenting large language models with chemistry tools, Nature Machine Intelligence 6 (2024) 525–535. doi:10.1038/s42256-024-00832-8. [28] D. Park, H. Moon, S. Ryu, A self-correcting multi-agent llm framework for language-based physics simulation and explanation, npj Artificial Intelligence 2 (2026) 10. doi:10.1038/s44387-025-00057-z. [29] K. Stechly, K. Valmeekam, S. Kambhampati, On the self-verification limitations of large language models on reasoning and planning tasks, proceedings of ICLR 2025. arXiv:2402.08115 (2025). [30] G. Dolcetti, E. Iotti, A dual perspective review on large language models and code verification, Frontiers in Computer Science 7 (2025) 1655469. doi:10.3389/fcomp.2025.1655469. [31] N. Tihanyi, T. Bisztray, M. A. Ferrag, R. Jain, L. C. Cordeiro, How secure is ai-generated code: A large-scale comparison of large language models, Empirical Software Engineering 30 (2) (2025) 47. doi:10.1007/s10664-024-10590-1. [32] M. N. Tran, M. F. Randolph, D. W. Airey, Study of seepage flow and sand plug loosening in installation of suction caissons in sand, in: Proceedings of the Fifteenth International Offshore and Polar Engineering Conference, International Society of Offshore and Polar Engineers, Seoul, Korea, 2005, p. 516–521. [33] L. Xie, S. Ma, T. Lin, The seepage and soil plug formation in suction caissons in sand using visual tests, Applied Sciences 10 (2) (2020) 566. doi:10.3390/app10020566. [34] K. Been, M. G. Jefferies, A state parameter for sands, Géotechnique 35 (2) (1985) 99–112. doi:10.1680/geot.1985.35.2.99. [35] J. Jaky, The coefficient of earth pressure at rest, Journal of the Society of Hungarian Architects and Engineers 78 (22) (1944) 355–358, in Hungarian. [36] P. W. Mayne, F. H. Kulhawy, k0 29035_ 28720–ocr relationships in soil, Journal of the Geotechnical Engineering Division 108 (6) (1982) 851–872. doi:10.1061/AJGEB6.0001306. [37] L. P. Suwal, R. Kuwano, Statically and dynamically measured poisson’s ratio of granular soils on triaxial laboratory specimens, Geotechnical Testing Journal 36 (4) (2013) 493–505. doi:10.1520/GTJ20120108. [38] X. S. Li, Y. F. Dafalias, Dilatancy for cohesionless soils, Géotechnique 50 (4) (2000) 449–460. doi:10.1680/geot.2000.50.4.449. [39] M. G. Jefferies, Nor-sand: a simple critical state model for sand, Géotechnique 43 (1) (1993) 91–103. doi:10.1680/geot.1993.43.1.91. [40] M. D. Bolton, The strength and dilatancy of sands, Géotechnique 36 (1) (1986) 65–78. doi:10.1680/geot.1986.36.1.65. [41] T. I. Tjelta, Geotechnical experience from the installation of the europipe jacket with bucket foundations, in: Offshore Technology Conference, 1995, p. OTC 7795. doi:10.4043/7795-MS. [42] A. K. Koteras, L. B. Ibsen, Medium-scale laboratory model of mono-bucket foundation for installation tests in sand, Canadian Geotechnical Journal 56 (8) (2019) 1142–1153. doi:10.1139/cgj-2018-0134. [43] J. Zhao, D. Li, Y. Zhang, Y. Zhang, Formation of soil plug in bioinspired scaled suction caissons in sand for offshore wind turbine, Energy 332 (2025) 137157. doi:10.1016/j.energy.2025.137157. [44] Y. Wang, X. Zhu, Q. Yang, Installation of suction caissons and formation of soil plug considering variation of permeability of sand, Chinese Journal of Geotechnical Engineering 41 (1) (2019) 184–190. doi:10.11779/CJGE201901021. [45] Q. Jiang, G. Karniadakis, AgenticSciML: collaborative multi-agent systems for emergent discovery in scientific machine learning, npj Artificial Intelligence 2 (2026) 57. doi:10.1038/s44387-026-00102-5. [46] A. E. Ghareeb, B. Chang, L. Mitchener, A. Yiu, C. J. Szostkiewicz, D. Shved, G. J. Gyimesi, J. M. Laurent, S. M. Wright, M. T. Razzak, A. D. White, S. C. Finnemann, M. M. Hinks, S. G. Rodriques, A multi-agent system for automating scientific discovery, Nature 655 (2026) 497–505. doi:10.1038/s41586-026-10652-y. [47] S. K. Thota, T. D. Cao, F. Vahedifard, Poisson’s ratio characteristic curve of unsaturated soils, Journal of Geotechnical and Geoenvironmental Engineering 147 (1) (2021) 04020149. doi:10.1061/(ASCE)GT.1943-5606.0002424. [48] S. I. Woo, C. K. Chung, Numerical simulation of cone penetration tests in sand ground using critical state mohr coulomb plasticity model, Journal of the Korean Geotechnical Society 35 (2) (2019) 37–51. doi:10.7843/kgs.2019.35.2.37. [49] Z. Zhu, W. Cheng, Parameter evaluation of exponential-form critical state line of a state-dependent sand constitutive model, Applied Sciences 10 (1) (2020) 328. doi:10.3390/app10010328. [50] W. Liu, M. Ghafghazi, Evaluation of state parameter interpretation methods using cpt calibration chamber data, Geotechnical Research 11 (4) (2024) 226–242. doi:10.1680/jgere.23.00072. [51] R. Kurniadi, A. Roy, S. H. Chow, M. J. Cassidy, Drained uplift capacity of deeply embedded horizontal circular plate anchors in sand, Ocean Engineering 331 (2025) 121259. doi:10.1016/j.oceaneng.2025.121259. Funding J.S. and Z.O. acknowledge support from the National Natural Science Foundation of China (grant nos. 52401333 and U24B20114). Author contributions J.S. conducted the research and wrote the original draft. J.S. and Y.L. prepared the figures. Y.L. revised the manuscript. Z.O. supervised the research and acquired funding. Ethics declaration Not applicable. Competing interests The authors declare no competing interests. Supplementary Information Human-guided physics-constrained AI agents construct an auditable model of soil-plug evolution This Supplementary Information first presents the model equations and numerical implementation of the soil-plug solver. It then reports the procedures and findings of the two workflow audits, followed by the checks applied to the final code. The final note and accompanying data files document the benchmark inputs, their sources, and the numerical results cited in the main text. Contents Supplementary Table 1. Contents of the Supplementary Information. Section Content Supplementary Note 1 Model formulation and numerical solution procedure. Supplementary Note 2 Design and findings of the two workflow audits. Supplementary Note 3 Verification of the final solver through module-level, regression, and analytical reference tests. Supplementary Note 4 Benchmark inputs, data sources, and comparison results. Supplementary Data 1 Numerical data underlying the benchmark summaries and Fig. 5 of the main text. Supplementary Data 2 Sources and uses of the global and case-specific model inputs. The benchmark evaluation reported in the main text uses the boundary rules, global settings, and study-specific inputs documented in Supplementary Tables 10–12. Supplementary Tables 13–16 report the three corresponding comparisons shown in Fig. 5. Supplementary Note 1. Mathematical Formulation and Numerical Implementation of the Soil-Plug Model Supplementary Note 1 defines how the soil-plug formulation selected in the main text calculates heave during caisson penetration and how this calculation is implemented numerically. The formulation begins with solid-volume conservation, which combines geometric displacement, seepage-induced changes in soil volume, and near-wall dilation. The following sections derive these contributions, assemble the final heave calculation, and then present the discrete solution process, numerical checks, and range of application. S1.1 Overall Formulation and Solution Target The soil-plug model calculates heave by combining three physical contributions. MG 29005_ 28999 represents the geometric displacement caused by skirt penetration, MS 29005_ 29011 represents the change in soil void ratio driven by seepage and effective-stress changes, and MD 29005_ 28996 represents near-wall soil dilation. At penetration depth z 29050, solid-volume conservation combines these contributions as follows: hea|e(z)=∫0z[αA1+e(ζ0,z)1+e0−1]dζ0, 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 12349 4946 _ 28720 29050 67482370 28939_ 28993 28721 8235 29029 67273472 28944_ 28720 24891 29050 84054785 28721 8235 29029_ 28720 8704 28721 84267779\, 29028 28944_ 28720 24891 (S1) where ζ0 28944_ 28720 labels a soil layer by its initial depth below the seabed, e(ζ0,z) 29029 67273472 28944_ 28720 24891 29050 84054785 is its void ratio after the MS 29005_ 29011 update and subsequent MD 29005_ 28996 correction, and e0 29029_ 28720 is the initial void ratio. The area ratio is αA=(Do/Di)2 28939_ 28993 12349 67273472 28996_ 29039 68408078 28996_ 29033 84054785 28722, where Do 28996_ 29039 and Di 28996_ 29033 are the outer and inner caisson diameters. Equation (S1) calculates plug heave from the void-ratio profile inside the caisson. MG 29005_ 28999 first establishes the geometric limit in which the void ratio remains unchanged. To determine the departure from that limit, MS 29005_ 29011 converts the prescribed suction-pressure history into a hydraulic gradient, an effective-stress path, and then a seepage-driven void-ratio update. After the MG 29005_ 28999–MS 29005_ 29011 calculation has converged, MD 29005_ 28996 adds the near-wall dilation mobilized during the penetration increment. The final void-ratio field is then returned to Eq. (S1) to calculate plug length and plug heave. Supplementary Table 2. Symbols, parameters, and state conventions used in the soil-plug formulation. Symbol or mark Definition Unit or type Geometry and coordinates Di,Do,ri 28996_ 29033 24891 28996_ 29039 24891 29042_ 29033 Di 28996_ 29033: inner caisson diameter; Do 28996_ 29039: outer caisson diameter; ri 29042_ 29033: inner caisson radius. m Ai,αA 28993_ 29033 24891 28939_ 28993 Ai 28993_ 29033: internal cross-sectional area; αA 28939_ 28993: outer-to-inner area ratio, (Do/Di)2 67273472 28996_ 29039 68408078 28996_ 29033 84054785 28722. Ai 28993_ 29033: m2; αA 28939_ 28993: dimensionless z,zn,zn+1 29050 24891 29050_ 29038 24891 29050_ 29038 8235 28721 z 29050: penetration depth; zn 29050_ 29038: depth at the start of the current step; zn+1 29050_ 29038 8235 28721: depth at the end of the current step. m Δz,z˙,Δt 28673 29050 24891 29050 24891 28673 29044 Δz 28673 29050: penetration increment; z˙ 29050: penetration rate; Δt 28673 29044: elapsed time for one step. Δz 28673 29050: m; z˙ 29050: m s-1; Δt 28673 29044: s hea|e,Lplug,H 29032_ 29032 29029 29025 69640972 29029 24891 29004_ 29040 29036 29045 29031 24891 29000 hea|e 29032_ 29032 29029 29025 69640972 29029: net plug heave; Lplug 29004_ 29040 29036 29045 29031: total plug length; H 29000: domain height measured from the plug free surface, H=Lplug=z+hea|e 29000 12349 29004_ 29040 29036 29045 29031 12349 29050 8235 29032_ 29032 29029 29025 69640972 29029. m y,x 29049 24891 29048 y 29049: fixed vertical coordinate measured downward from the initial seabed; x 29048: vertical coordinate measured downward from the current plug free surface. m Hydraulic parameters and flow rates Δu 28673 29045 Prescribed suction-pressure difference. kPa γw,γ′ 28941_ 29047 24891 28941 560 γw 28941_ 29047: unit weight of water; γ′ 28941 560: buoyant unit weight of sand. kN m-3 ϕ 28958 Equivalent head loss. m i,itip 29033 24891 29033_ 29044 29033 29040 i 29033: hydraulic gradient, dϕ/dx 29028 28958 68408078 29028 29048; itip 29033_ 29044 29033 29040: hydraulic gradient at the caisson tip. dimensionless k|,kr 29035_ 69640972 24891 29035_ 29042 k| 29035_ 69640972: equivalent vertical permeability coefficient; kr 29035_ 29042: equivalent radial permeability coefficient. m s-1 ki,kout 29035_ 29033 24891 29035_ 29039 29045 29044 ki 29035_ 29033: internal-soil permeability; kout 29035_ 29039 29045 29044: external-soil permeability. m s-1 Re∗,Re,∞,Rb 29010_ 29029 8707 24891 29010_ 29029 24891 561 24891 29010_ 29026 Re∗ 29010_ 29029 8707: adopted equivalent outer-domain radius; Re,∞ 29010_ 29029 24891 561: modal radius estimate for an unbounded outer domain; Rb 29010_ 29026: source-based finite limiting radius. m ℓs,λ1 352_ 29043 24891 28949_ 28721 ℓs 352_ 29043: characteristic seepage length; λ1 28949_ 28721: first vertical modal wavenumber. ℓs 352_ 29043: m; λ1 28949_ 28721: m-1 K0(⋅),K1(⋅) 29003_ 28720 67273472 8705 84054785 24891 29003_ 28721 67273472 8705 84054785 K0(⋅) 29003_ 28720 67273472 8705 84054785: zeroth-order modified Bessel function of the second kind; K1(⋅) 29003_ 28721 67273472 8705 84054785: first-order modified Bessel function of the second kind. functions Qx,qr 29009_ 29048 24891 29041_ 29042 Qx 29009_ 29048: signed axial Darcy discharge; qr 29041_ 29042: radial seepage rate per unit plug length. Qx 29009_ 29048: m3 s-1; qr 29041_ 29042: m2 s-1 Qin,top,Qtip 29009_ 29033 29038 24891 29044 29039 29040 24891 29009_ 29044 29033 29040 Qin,top 29009_ 29033 29038 24891 29044 29039 29040: upward flow rate across the plug top; Qtip 29009_ 29044 29033 29040: upward flow rate across the caisson-tip section. m3 s-1 Qgeom,Qpump 29009_ 29031 29029 29039 29037 24891 29009_ 29040 29045 29037 29040 Qgeom 29009_ 29031 29029 29039 29037: geometric volume-change rate; Qpump 29009_ 29040 29045 29037 29040: equivalent pump flow. m3 s-1 Stress and soil-state parameters σ|′,σh′,σrc′ 28955 560_ 69640972 24891 28955 560_ 29032 24891 28955 560_ 29042 29027 σ|′ 28955 560_ 69640972: vertical effective stress; σh′ 28955 560_ 29032: lateral effective stress; σrc′ 28955 560_ 29042 29027: interfacial reference confining stress. kPa σ^′|,σ′|,reg 28955 560_ 69640972 24891 28955 560_ 69640972 24891 29042 29029 29031 σ^|′ 28955 560_ 69640972: regularized vertical effective stress in the refreshed final output; σ′|,reg 28955 560_ 69640972 24891 29042 29029 29031: regularized vertical effective stress in trial or history states. kPa p′,p^′,q 29040 560 24891 29040 560 24891 29041 p′ 29040 560: mean effective stress; p^′ 29040 560: regularized mean effective stress; q 29041: deviatoric stress magnitude. kPa η 28945 Stress ratio, q/p^′ 29041 68408078 29040 560. dimensionless ν,K0,Kp 28951 24891 29003_ 28720 24891 29003_ 29040 ν 28951: Poisson’s ratio; K0 29003_ 28720: at-rest earth-pressure coefficient; Kp 29003_ 29040: passive stress-ratio bound. dimensionless cp(K0) 29027_ 29040 67273472 29003_ 28720 84054785 Coefficient converting vertical effective stress to the stored mean-stress reference, (1+2K0)/3 67273472 28721 8235 28722 29003_ 28720 84054785 68408078 28723. dimensionless pmin′ 29040 560_ Positive numerical stress floor. kPa ϕmodel 28958_ 29037 29039 29028 29029 29036 Case-specific friction-angle input. degrees e,e0 29029 24891 29029_ 28720 e 29029: current void ratio; e0 29029_ 28720: initial void ratio. dimensionless emin,emax 29029_ 24891 29029_ emin 29029_ : source-based lower void-ratio bound; emax 29029_ : source-based upper void-ratio bound. dimensionless ecs,Ψ 29029_ 29027 29043 24891 28681 ecs 29029_ 29027 29043: critical-state void ratio; Ψ 28681: state parameter, e−ecs 29029 8704 29029_ 29027 29043. dimensionless eUNKNOWN,λc,κc 29029_ 28672 24891 28949_ 29027 24891 28948_ 29027 eUNKNOWN 29029_ 28672: critical-state-line intercept; λc 28949_ 29027: critical-state-line shape coefficient; κc 28948_ 29027: critical-state-line shape exponent. dimensionless pa 29040_ 29025 Reference stress in the critical-state relation. kPa Mtc,Mte,Mpath 29005_ 29044 29027 24891 29005_ 29044 29029 24891 29005_ 29040 29025 29044 29032 Mtc 29005_ 29044 29027: triaxial-compression critical stress ratio; Mte 29005_ 29044 29029: triaxial-extension critical stress ratio; Mpath 29005_ 29040 29025 29044 29032: critical stress ratio adopted for the current path. dimensionless ηact,η0 28945_ 29025 29027 29044 24891 28945_ 28720 ηact 28945_ 29025 29027 29044: stress ratio at node activation; η0 28945_ 28720: stored activation reference stress ratio. dimensionless μ,μ¯,Δμ¯ 28950 24891 28950 24891 28673 28950 μ 28950: current path-mobilization degree; μ¯ 28950: maximum mobilization reached in accepted steps; Δμ¯ 28673 28950: newly mobilized increment. dimensionless ϵUNKNOWN 28943_ 28945 Numerical denominator floor in the mobilization calculation. dimensionless κs 28948_ 29043 Swelling index. dimensionless Δereb,ereb 28673 29029_ 29042 29029 29026 24891 29029 29042 29029 29026 Δereb 28673 29029_ 29042 29029 29026: void-ratio increment due to rebound; ereb 29029 29042 29029 29026: void ratio after rebound. dimensionless Creb,Δerel 28995 29042 29029 29026 24891 28673 29029_ 29042 29029 29036 Creb 28995 29042 29029 29026: remaining dense-state capacity for void-ratio increase after rebound; Δerel 28673 29029_ 29042 29029 29036: void-ratio increment released by path mobilization. dimensionless Near-wall dilation parameters Dr,Q 28996_ 29042 24891 29009 Dr 28996_ 29042: relative density; Q 29009: parameter in the Bolton-type dilation relation. dimensionless c ,pa,B 29027_ 28960 24891 29040_ 29025 24891 28994 c 29027_ 28960: coefficient converting relative dilatancy index to dilatancy angle; pa,B 29040_ 29025 24891 28994: reference pressure in the dilation relation. c 29027_ 28960: degrees per unit IR 29001_ 29010; pa,B 29040_ 29025 24891 28994: kPa IR,ψd 29001_ 29010 24891 28960_ 29028 IR 29001_ 29010: Bolton-type relative dilatancy index; ψd 28960_ 29028: local dilatancy angle. IR 29001_ 29010: dimensionless; ψd 28960_ 29028: rad δy,δcum 28942_ 29049 24891 28942_ 29027 29045 29037 δy 28942_ 29049: characteristic displacement; δcum 28942_ 29027 29045 29037: cumulative interfacial displacement recorded in the dilation history. m Δδpot,Δδmob 28673 28942_ 29040 29039 29044 24891 28673 28942_ 29037 29039 29026 Δδpot 28673 28942_ 29040 29039 29044: potential interfacial displacement increment; Δδmob 28673 28942_ 29037 29039 29026: increment recorded when a positive dilation increment is applied. m Δε¯|,maxD,Δε¯|,rawD 28673 28962_ 69640972 24891 28996 24891 28673 28962_ 69640972 24891 29042 29025 29047 28996 Δε¯|,maxD 28673 28962_ 69640972 24891 28996: maximum cross-sectionally averaged dilation strain increment; Δε¯|,rawD 28673 28962_ 69640972 24891 29042 29025 29047 28996: increment before capacity truncation. dimensionless Δε¯|,remD,Δε¯|D 28673 28962_ 69640972 24891 29042 29029 29037 28996 24891 28673 28962_ 69640972 28996 Δε¯|,remD 28673 28962_ 69640972 24891 29042 29029 29037 28996: remaining cross-sectionally averaged dilation strain increment; Δε¯|D 28673 28962_ 69640972 28996: applied increment. dimensionless εcumD,ω 28962_ 29027 29045 29037 28996 24891 28961 εcumD 28962_ 29027 29045 29037 28996: accumulated dilation strain amplitude; ω 28961: dilation-capacity decay factor. dimensionless ΔeD 28673 29029 28996 Void-ratio increment due to dilation. dimensionless Discretization, iteration, and state notation ζ0,j,Δζ0,j 28944_ 28720 24891 29034 24891 28673 28944_ 28720 24891 29034 ζ0,j 28944_ 28720 24891 29034: initial depth of material node j 29034; Δζ0,j 28673 28944_ 28720 24891 29034: trapezoidal quadrature weight assigned to that node. m j,N 29034 24891 29006 j 29034: material-node index; N 29006: total number of material nodes. integers n,k,kcon| 29038 24891 29035 24891 29035_ 29027 29039 29038 69640972 n 29038: penetration-step index; k 29035: fixed-point iteration index; kcon| 29035_ 29027 29039 29038 69640972: accepted converged iteration index. integers hest,hcand 29032_ 29029 29043 29044 24891 29032_ 29027 29025 29038 29028 hest 29032_ 29029 29043 29044: plug-heave estimate entering an iteration; hcand 29032_ 29027 29025 29038 29028: candidate plug heave calculated in that iteration. m rh 29042_ 29032 Absolute plug-heave residual used to test fixed-point convergence. m ℐjn+1,n+1 29001_ 29034 29038 8235 28721 24891 28993_ 29038 8235 28721 ℐjn+1 29001_ 29034 29038 8235 28721: node-activation flag for the current step; n+1 28993_ 29038 8235 28721: active-node set for the current step. ℐjn+1 29001_ 29034 29038 8235 28721: 0 or 1; n+1 28993_ 29038 8235 28721: set j,n+1∗ 29011_ 29034 24891 29038 8235 28721 8707 Synchronized state of node j 29034 after coupled convergence and before dilation correction. state vector ℱn+1 28998_ 29038 8235 28721 Fixed-point map from trial void ratio and plug heave to their candidate values. map tr,H 29044 29042 24891 29000 tr 29044 29042: trial-state mark; H 29000: stored-history mark. superscripts D,∗ 28996 24891 8707 D 28996: dilation-stage mark; ∗ 8707: synchronized-state mark after coupled convergence and before dilation correction. superscripts out,con| 29039 29045 29044 24891 29027 29039 29038 69640972 out 29039 29045 29044: refreshed final-output mark; con| 29027 29039 29038 69640972: converged-iteration mark. superscripts or labels proj,sat,clip 29040 29042 29039 29034 24891 29043 29025 29044 24891 29027 29036 29033 29040 proj 29040 29042 29039 29034: projection to stress bounds; sat 29043 29025 29044: restriction to [0,1] 67482370 28720 24891 28721 84267779; clip 29027 29036 29033 29040: restriction to an admissible interval. operators Model components MG,MS,MD 29005_ 28999 24891 29005_ 29011 24891 29005_ 28996 MG 29005_ 28999: skirt-penetration-induced geometric displacement; MS 29005_ 29011: seepage-driven soil void-ratio evolution; MD 29005_ 28996: near-wall soil dilation. model components S1.2 Coordinate System and Geometric Baseline Let z 29050 denote the caisson penetration depth and hea|e(z) 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 the net plug heave. The fixed coordinate y 29049 is measured downward from the initial mudline. The coordinate attached to the moving plug free surface is x=y+hea|e(z),H=z+hea|e(z), 29048 12349 29049 8235 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 24891 29000 12349 29050 8235 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 24891 (S2) where x 29048 is positive downward, x=0 29048 12349 28720 denotes the current plug free surface, x=H 29048 12349 29000 denotes the caisson tip, and H=Lplug(z) 29000 12349 29004_ 29040 29036 29045 29031 67273472 29050 84054785 is the current total plug length. When the void ratio remains at its initial value, e(ζ0,z)=e0 29029 67273472 28944_ 28720 24891 29050 84054785 12349 29029_ 28720, Eq. (S1) reduces to hea|e(z)=(αA−1)z. 29032_ 29032 29029 29025 69640972 29029 67273472 29050 84054785 12349 67273472 28939_ 28993 8704 28721 84054785 29050 314 (S3) This relation defines the MG 29005_ 28999 geometric baseline. The MS 29005_ 29011 and MD 29005_ 28996 components move the model away from this baseline by changing the void-ratio field used in the solid-volume-conservation relation. S1.3 Seepage-Driven Void-Ratio Evolution (MS 29005_ 29011) Prescribed suction first creates a hydraulic head difference across the soil plug. This difference determines the hydraulic gradient, which changes the effective stresses carried by the soil skeleton. The resulting stress path then determines the change in void ratio. MS 29005_ 29011 therefore follows the sequence Δu→i→(σ|′,σh′,p′,q)→e∗ 28673 29045 12833 29033 12833 67273472 28955 560_ 69640972 24891 28955 560_ 29032 24891 29040 560 24891 29041 84054785 12833 29029 8707, where e∗ 29029 8707 is the void ratio accepted after the coupled MG 29005_ 28999–MS 29005_ 29011 calculation converges and then passed to MD 29005_ 28996. S1.3.1 From Prescribed Suction to the Hydraulic Field The first step converts the prescribed suction-pressure difference into the hydraulic gradient along the plug. The hydraulic calculation resolves vertical Darcy flow along x 29048 and represents seepage through the surrounding soil as an equivalent radial discharge per unit plug length. The two flow paths form the one-dimensional axial equation used below. Let ϕ(x) 28958 67273472 29048 84054785 denote the equivalent head loss, expressed as water head, and let Δu(z)>0 28673 29045 67273472 29050 84054785 12606 28720 denote the prescribed suction-pressure difference. With γw 28941_ 29047 denoting the unit weight of water, the boundary conditions are ϕ(0)=0,ϕ(H)=Δu(z) w. 28958 67273472 28720 84054785 12349 28720 24891 28958 67273472 29000 84054785 12349 28673 29045 67273472 29050 84054785 28941_ 29047 314 (S4) The signed axial discharge Qx 29009_ 29048 uses the equivalent vertical permeability k| 29035_ 69640972 and the internal cross-sectional area Ai 28993_ 29033: Qx=−k|Aidϕdx 29009_ 29048 12349 8704 29035_ 69640972 28993_ 29033 29028 28958 29028 29048 (S5) The radial seepage rate per unit plug length, qr 29041_ 29042, uses the equivalent radial permeability kr 29035_ 29042, inner radius ri 29042_ 29033, and adopted equivalent outer-domain radius Re∗ 29010_ 29029 8707. In this equivalent one-dimensional plug domain, qr>0 29041_ 29042 12606 28720 denotes radial discharge leaving the plug domain: qr=2πkrϕln(Re∗/ri). 29041_ 29042 12349 28722 28953 29035_ 29042 28958 (R_e^*/r_i) 314 (S6) Conservation of fluid mass gives dQx/dx+qr=0 29028 29009_ 29048 68408078 29028 29048 8235 29041_ 29042 12349 28720. With Ai=πri2 28993_ 29033 12349 28953 29042_ 29033 28722, this balance becomes d2ϕdx2−2krk|ri2ln(Re∗/ri)ϕ=0. 29028 28722 28958 29028 29048 28722 8704 28722 29035_ 29042 29035_ 69640972 29042_ 29033 28722 (R_e^*/r_i) 28958 12349 28720 314 (S7) The radius Re∗ 29010_ 29029 8707 therefore controls the radial drainage resistance represented in the one-dimensional model. Defining ℓs 352_ 29043 as the characteristic seepage length gives the equivalent one-dimensional finite-domain kernel d2ϕdx2−ϕℓs2=0, 29028 28722 28958 29028 29048 28722 8704 28958 352_ 29043 28722 12349 28720 24891 (S8) with ri=Di2,ℓs=ri[k|2krln(Re∗ri)]1/2, 29042_ 29033 12349 28996_ 29033 28722 24891 352_ 29043 12349 29042_ 29033 67482370 29035_ 69640972 28722 29035_ 29042 67273472 29010_ 29029 8707 29042_ 29033 84054785 84267779 28721 68408078 28722 24891 (S9) where Di 28996_ 29033 is the inner diameter. For each benchmark case, we set the equivalent ratio kr/k| 29035_ 29042 68408078 29035_ 69640972 in Eq. (S9) equal to the adopted internal-to-external permeability ratio ki/kout 29035_ 29033 68408078 29035_ 29039 29045 29044 as a model assumption. The unbounded modal estimate Re,∞ 29010_ 29029 24891 561 is defined by matching the logarithmic radial drainage resistance of the primary leakage mode, ln(Re,∞ri)=K0(λ1ri)UNKNOWN1riK1(UNKNOWN1ri),λ1=πH, 67273472 29010_ 29029 24891 561 29042_ 29033 84054785 12349 29003_ 28720 67273472 28949_ 28721 29042_ 29033 84054785 28949_ 28721 29042_ 29033 29003_ 28721 67273472 28949_ 28721 29042_ 29033 84054785 24891 28949_ 28721 12349 28953 29000 24891 (S10) where λ1 28949_ 28721 is the first vertical modal wavenumber and K0(⋅) 29003_ 28720 67273472 8705 84054785 and K1(⋅) 29003_ 28721 67273472 8705 84054785 are modified Bessel functions of the second kind. When a finite external boundary or source-based limiting radius Rb 29010_ 29026 is available, the adopted radius is Re∗=min(Rb,Re,∞) 29010_ 29029 8707 12349 67273472 29010_ 29026 24891 29010_ 29029 24891 561 84054785, with Re∗>ri 29010_ 29029 8707 12606 29042_ 29033. If no finite limiting radius is available, Re∗=Re,∞ 29010_ 29029 8707 12349 29010_ 29029 24891 561. The limiting radius is specified from source boundary information or from the benchmark parameter rule. Solving Eq. (S8) with Eq. (S4) gives ϕ(x)=Δu(z) wsinh(x/ℓs)sinh(H/ℓs). 28958 67273472 29048 84054785 12349 28673 29045 67273472 29050 84054785 28941_ 29047 (x/ _s) (H/ _s) 314 (S11) In the deep-plug limit, and away from the immediate free-surface boundary, Eq. (S11) approaches ϕ(x)≈Δu(z) wexp(x−Hℓs). 28958 67273472 29048 84054785 12825 28673 29045 67273472 29050 84054785 28941_ 29047 67273472 29048 8704 29000 352_ 29043 84054785 314 (S12) The hydraulic gradient is i=dϕ/dx 29033 12349 29028 28958 68408078 29028 29048; positive i 29033 represents the upward seepage contribution that reduces vertical effective stress: i(x,z)=Δu(z) wℓsinh(H/ℓs)cosh(x/ℓs). 29033 67273472 29048 24891 29050 84054785 12349 28673 29045 67273472 29050 84054785 28941_ 29047 352_ 29043 (H/ _s) (x/ _s) 314 (S13) The kernel satisfies the imposed pressure-difference constraint: γw∫0Hi(x,z)x=Δu(z). 28941_ 29047 4946 _ 28720 29000 29033 67273472 29048 24891 29050 84054785\, 29028 29048 12349 28673 29045 67273472 29050 84054785 314 (S14) The hydraulic gradient in Eq. (S13) supplies the seepage-force term required to calculate the vertical effective stress. S1.3.2 From Hydraulic Gradient to Effective Stress Upward seepage reduces the vertical effective stress carried by the soil skeleton. Combining the submerged soil weight with the seepage-force term gives the vertical effective stress at a position x 29048 below the current plug free surface: σ|′(x,z)=∫0x[γ′−γwi(s,z)]s, 28955 560_ 69640972 67273472 29048 24891 29050 84054785 12349 4946 _ 28720 29048 67482370 28941 560 8704 28941_ 29047 29033 67273472 29043 24891 29050 84054785 84267779 29028 29043 24891 (S15) where γ′ 28941 560 is the buoyant unit weight of sand. The symbol s 29043 is a temporary variable used to sum the soil-weight and seepage-force contributions from the plug surface to position x 29048. Substituting Eq. (S13) gives σ|′(x,z)=γ′x−Δu(z)sinh(x/ℓs)sinh(H/ℓs). 28955 560_ 69640972 67273472 29048 24891 29050 84054785 12349 28941 560 29048 8704 28673 29045 67273472 29050 84054785 (x/ _s) (H/ _s) 314 (S16) The free-surface boundary condition is σ|′(0,z)=0. 28955 560_ 69640972 67273472 28720 24891 29050 84054785 12349 28720 314 (S17) The hydraulic gradient at the caisson tip is obtained by evaluating Eq. (S13) at x=H 29048 12349 29000: itip(z)=i(H,z)=Δu(z) wℓscoth(H/ℓs). 29033_ 29044 29033 29040 67273472 29050 84054785 12349 29033 67273472 29000 24891 29050 84054785 12349 28673 29045 67273472 29050 84054785 28941_ 29047 352_ 29043 67273472 29000 68408078 352_ 29043 84054785 314 (S18) The resulting vertical stress field is then used to calculate the lateral, mean, and deviatoric stresses. The corresponding stress path controls the void-ratio update. S1.3.3 From Effective-Stress Path to Void-Ratio Change At each penetration depth, we use the effective stresses to calculate changes in void ratio and obtain plug heave from Eq. (S1). The effective stresses depend on seepage through the plug, which is itself affected by the plug height. Updating the void ratio thus changes the height used to calculate seepage and effective stresses. We account for this dependence by repeatedly recalculating seepage and stresses with the updated plug geometry, then updating the void ratio and heave. This procedure is implemented as a fixed-point iteration. Within each iteration, the void-ratio update combines two components. Changes in mean effective stress produce a rebound component, with unloading increasing the void ratio. A second component allows a further increase toward the critical-state void ratio as the stress ratio approaches its limiting value. For penetration step n+1 29038 8235 28721, k 29035 identifies the current iteration. The trial geometry defines the total plug length Hn+1(k) 29000_ 29038 8235 28721 67273472 29035 84054785 and the free-surface-attached coordinate xj,n+1(k) 29048_ 29034 24891 29038 8235 28721 67273472 29035 84054785 of node j 29034. The corresponding vertical effective stress is σ|,j,n+1′(k)=γ′xj,n+1(k)−Δu(zn+1)sinh(xj,n+1(k)/ℓs(k))sinh(Hn+1(k)/ℓs(k)). 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12349 28941 560 29048_ 29034 24891 29038 8235 28721 67273472 29035 84054785 8704 28673 29045 67273472 29050_ 29038 8235 28721 84054785 67273472 29048_ 29034 24891 29038 8235 28721 67273472 29035 84054785 68408078 352_ 29043 67273472 29035 84054785 84054785 67273472 29000_ 29038 8235 28721 67273472 29035 84054785 68408078 352_ 29043 67273472 29035 84054785 84054785 314 (S19) Nodes with σ|,j,n+1′(k)≤0 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12820 28720 are flagged as near-critical seepage locations. To keep the logarithmic and stress-ratio calculations well defined, the current trial stress and the stress stored from the preceding accepted step are bounded below by pmin′ 29040 560_ : σ|,j,n+1,reg′(k) 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 =max[σ|,j,n+1′(k),p′min], 12349 67482370 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 24891 29040 560_ 84267779 24891 (S20) σ|,j,n,reg′H 28955_ 69640972 24891 29034 24891 29038 24891 29042 29029 29031 560 29000 =max[σ|,j,n′H,p′min], 12349 67482370 28955_ 69640972 24891 29034 24891 29038 560 29000 24891 29040 560_ 84267779 24891 where pmin′ 29040 560_ is a positive numerical stress floor and the superscript H 29000 denotes a history value stored from the preceding accepted step. The lateral effective stress is first estimated from the change in regularized vertical stress relative to the preceding accepted state: σh,j,n+1′tr,(k)=σh,j,n′H+ν1−UNKNOWN(σ|,j,n+1,reg′(k)−σ|,j,n,reg′H), 28955_ 29032 24891 29034 24891 29038 8235 28721 560 29044 29042 24891 67273472 29035 84054785 12349 28955_ 29032 24891 29034 24891 29038 560 29000 8235 28951 28721 8704 28951 67273472 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 8704 28955_ 69640972 24891 29034 24891 29038 24891 29042 29029 29031 560 29000 84054785 24891 (S21) where ν 28951 is Poisson’s ratio. The trial value is then limited between the numerical stress floor and the passive stress bound: σh,j,n+1′(k)=proj[σh,j,n+1′tr,(k);pmin′,Kpσ|,j,n+1,reg′(k)]. 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12349 29040 29042 29039 29034 67482370 28955_ 29032 24891 29034 24891 29038 8235 28721 560 29044 29042 24891 67273472 29035 84054785 24635 29040 560_ 24891 29003_ 29040 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 84267779 314 (S22) The projection operator and passive stress-ratio bound are defined as proj[a;l,u]=min[max(a,l),u],Kp=1+sinϕmodel1−sinUNKNOWNmodel, 29040 29042 29039 29034 67482370 29025 24635 29036 24891 29045 84267779 12349 67482370 67273472 29025 24891 29036 84054785 24891 29045 84267779 24891 29003_ 29040 12349 28721 8235 28958_ 29037 29039 29028 29029 29036 28721 8704 28958_ 29037 29039 29028 29029 29036 24891 (S23) where a 29025 is the value being projected and l 29036 and u 29045 are its lower and upper bounds. The quantity Kp 29003_ 29040 is the passive stress-ratio bound, and ϕmodel 28958_ 29037 29039 29028 29029 29036 is the case-specific friction-angle input used by the equivalent one-dimensional solver. The vertical and lateral effective stresses are used to calculate the mean effective stress and the stress ratio. Of these, the mean effective stress is used to determine the rebound response and the critical-state void ratio, while the stress ratio is used to quantify progress toward the limiting stress state. The mean effective stress p′ 29040 560 and the corresponding regularized value p^′ 29040 560 are pj,n+1′(k)=σ|,j,n+1,reg′(k)+2σh,j,n+1′(k)3,p^j,n+1′(k)=max[pj,n+1′(k),pmin′]. 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12349 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 8235 28722 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 28723 24891 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12349 67482370 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 24891 29040 560_ 84267779 314 (S24) The deviatoric stress magnitude q 29041 and stress ratio η=q/p^′ 28945 12349 29041 68408078 29040 560 are qj,n+1(k)=|σ|,j,n+1,reg′(k)−σh,j,n+1′(k)|,ηj,n+1(k)=qj,n+1(k)p^j,n+1′(k). 29041_ 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 69640972 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 8704 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 69640972 24891 28945_ 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 29041_ 29034 24891 29038 8235 28721 67273472 29035 84054785 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 314 (S25) The critical-state void ratio provides the reference for the additional void-ratio increase. Its value at the current mean effective stress is ecs,j,n+1(k)=eUNKNOWN−λc(p^j,n+1′(k)pa)UNKNOWNc, 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 29029_ 28672 8704 28949_ 29027 67273472 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 29040_ 29025 84054785 28948_ 29027 24891 (S26) where eUNKNOWN 29029_ 28672, λc 28949_ 29027, and κc 28948_ 29027 are the intercept and shape parameters of the adopted critical-state line, and pa 29040_ 29025 is its reference stress. To measure progress toward the limiting stress state, the current stress ratio is compared with its value when the soil first enters the caisson and with the applicable stress-ratio limit. The entry value is stored for each node as η0,j=ηj,act. 28945_ 28720 24891 29034 12349 28945_ 29034 24891 29025 29027 29044 314 (S27) The friction-angle input defines separate stress-ratio limits for triaxial compression and extension: Mtc=6sinϕmodel3−sinUNKNOWNmodel,Mte=6sinϕmodel3+sinUNKNOWNmodel. 29005_ 29044 29027 12349 28726 28958_ 29037 29039 29028 29029 29036 28723 8704 28958_ 29037 29039 29028 29029 29036 24891 29005_ 29044 29029 12349 28726 28958_ 29037 29039 29028 29029 29036 28723 8235 28958_ 29037 29039 29028 29029 29036 314 (S28) The applicable limit depends on whether the vertical or lateral effective stress is larger: Mpath,j,n+1(k)=Mtc,UNKNOWN|,j,n+1,reg′(k)≥UNKNOWNh,j,n+1′(k),Mte,UNKNOWNh,j,n+1′(k)>UNKNOWN|,j,n+1,reg′(k). 29005_ 29040 29025 29044 29032 24891 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 cases 29005_ 29044 29027 24891& 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 12821 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 24891\\[4.0pt] 29005_ 29044 29029 24891& 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035 84054785 12606 28955_ 69640972 24891 29034 24891 29038 8235 28721 24891 29042 29029 29031 560 67273472 29035 84054785 314 cases (S29) The mobilization degree μ 28950 quantifies the increase in stress ratio from the entry value η0,j 28945_ 28720 24891 29034 toward the selected limit. The operator sat[a,0,1]=min[max(a,0),1] 29043 29025 29044 67482370 29025 24891 28720 24891 28721 84267779 12349 67482370 67273472 29025 24891 28720 84054785 24891 28721 84267779 bounds its argument a 29025 between 0 and 1: μj,n+1(k)=sat[ηj,n+1(k)−η0,jmax(Mpath,j,n+1(k)−UNKNOWN0,j,UNKNOWNUNKNOWN),0,1], 28950_ 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 29043 29025 29044 67482370 28945_ 29034 24891 29038 8235 28721 67273472 29035 84054785 8704 28945_ 28720 24891 29034 67273472 29005_ 29040 29025 29044 29032 24891 29034 24891 29038 8235 28721 67273472 29035 84054785 8704 28945_ 28720 24891 29034 24891 28943_ 28945 84054785 24891 28720 24891 28721 84267779 24891 (S30) where ϵUNKNOWN 28943_ 28945 is a positive lower bound used to keep the denominator nonzero. Values of μ=0 28950 12349 28720, 0<μ<1 28720 12604 28950 12604 28721, and μ=1 28950 12349 28721 represent no, partial, and full mobilization, respectively. The rebound contribution is then calculated from the change in mean effective stress relative to the preceding accepted step. It uses the swelling index κs 28948_ 29043 and the positive mean-stress reference pj,n′H 29040_ 29034 24891 29038 560 29000 stored at that step: Δereb,j(k)=−κsln(p^j,n+1′(k)pj,n′H),ej,n+1reb,(k)=ej,n+Δereb,j(k). 28673 29029_ 29042 29029 29026 24891 29034 67273472 29035 84054785 12349 8704 28948_ 29043 67273472 29040_ 29034 24891 29038 8235 28721 560 67273472 29035 84054785 29040_ 29034 24891 29038 560 29000 84054785 24891 29029_ 29034 24891 29038 8235 28721 29042 29029 29026 24891 67273472 29035 84054785 12349 29029_ 29034 24891 29038 8235 28673 29029_ 29042 29029 29026 24891 29034 67273472 29035 84054785 314 (S31) After the rebound update, the remaining capacity for an additional void-ratio increase is the positive gap between the critical-state void ratio and the value after rebound: Cj,n+1reb,(k)=max[ecs,j,n+1(k)−ej,n+1reb,(k),0]. 28995_ 29034 24891 29038 8235 28721 29042 29029 29026 24891 67273472 29035 84054785 12349 67482370 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 67273472 29035 84054785 8704 29029_ 29034 24891 29038 8235 28721 29042 29029 29026 24891 67273472 29035 84054785 24891 28720 84267779 314 (S32) To count each mobilization increase once, the current degree is compared with the maximum saved through the preceding accepted step. The positive difference gives the new mobilization increment: Δμ¯j,n+1(k)=max[μj,n+1(k)−μ¯j,n,0], 28673 28950_ 29034 24891 29038 8235 28721 67273472 29035 84054785 12349 67482370 28950_ 29034 24891 29038 8235 28721 67273472 29035 84054785 8704 28950_ 29034 24891 29038 24891 28720 84267779 24891 (S33) where μ¯j,n 28950_ 29034 24891 29038 is the maximum mobilization degree committed through step n 29038. Multiplying this new mobilization increment by the remaining capacity gives the additional void-ratio increase: Δerel,j(k)=Δμ¯j,n+1(k)Cj,n+1reb,(k). 28673 29029_ 29042 29029 29036 24891 29034 67273472 29035 84054785 12349 28673 28950_ 29034 24891 29038 8235 28721 67273472 29035 84054785 28995_ 29034 24891 29038 8235 28721 29042 29029 29026 24891 67273472 29035 84054785 314 (S34) Adding this increase to the void ratio after rebound gives the trial void ratio. The operator clip[a,l,u]=min[max(a,l),u] 29027 29036 29033 29040 67482370 29025 24891 29036 24891 29045 84267779 12349 67482370 67273472 29025 24891 29036 84054785 24891 29045 84267779 limits this value to the source-based admissible bounds emin 29029_ and emax 29029_ : ej,n+1∗,(k)=clip[ej,n+1reb,(k)+Δerel,j(k),emin,emax]. 29029_ 29034 24891 29038 8235 28721 8707 24891 67273472 29035 84054785 12349 29027 29036 29033 29040 67482370 29029_ 29034 24891 29038 8235 28721 29042 29029 29026 24891 67273472 29035 84054785 8235 28673 29029_ 29042 29029 29036 24891 29034 67273472 29035 84054785 24891 29029_ 24891 29029_ 84267779 314 (S35) The superscript ∗ 8707 denotes the coupled-loop state before the MD 29005_ 28996 correction. The state parameter Ψ=e−ecs 28681 12349 29029 8704 29029_ 29027 29043 measures the difference between this updated void ratio and its critical-state value: Ψj,n+1∗,(k)=ej,n+1∗,(k)−ecs,j,n+1(k). 28681_ 29034 24891 29038 8235 28721 8707 24891 67273472 29035 84054785 12349 29029_ 29034 24891 29038 8235 28721 8707 24891 67273472 29035 84054785 8704 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 67273472 29035 84054785 314 (S36) S1.4 Near-Wall Soil Dilation (MD 29005_ 28996) MD 29005_ 28996 calculates the additional void-ratio increment caused by near-wall soil dilation after the MG 29005_ 28999–MS 29005_ 29011 calculation has converged. To obtain this increment, MD 29005_ 28996 needs the local dilation tendency, the interface displacement mobilized during the current penetration step, and the dilation capacity remaining in the current soil state. These quantities jointly determine the void-ratio increment ΔejD 28673 29029_ 29034 28996, which is added to the converged void ratio. S1.4.1 Converged Input State and Near-Wall Confinement To calculate the local dilation tendency, MD 29005_ 28996 establishes the normal confinement acting near the caisson wall. For each active soil layer, the calculation begins from the void ratio and effective stresses returned by the converged MG 29005_ 28999–MS 29005_ 29011 calculation: ej∗,σ|,j′∗,σh,j′∗, \ 29029_ 29034 8707 24891\, 28955_ 69640972 24891 29034 560 8707 24891\, 28955_ 29032 24891 29034 560 8707 \ 24891 (S37) where the superscript ∗ 8707 marks the synchronized state accepted at MG 29005_ 28999–MS 29005_ 29011 convergence. The lateral effective stress provides the wall confinement used in the dilation relation. The reference confinement is bounded below by pmin′ 29040 560_ as a numerical safeguard: σrc,j′=max(σh,j′∗,pmin′). 28955 560_ 29042 29027 24891 29034 12349 67273472 28955_ 29032 24891 29034 560 8707 24891 29040 560_ 84054785 314 (S38) S1.4.2 Dilation Tendency and Potential Interface Displacement With the wall confinement established, MD 29005_ 28996 calculates whether the soil tends to dilate and the strength of that tendency. The adopted Bolton-type relation expresses this tendency as the relative dilatancy index IR,j 29001_ 29010 24891 29034 [40]. The relation evaluates the index from the relative density Dr 28996_ 29042, material parameter Q 29009, reference pressure pa,B 29040_ 29025 24891 28994, and wall confinement: IR,j=max[Dr(Q−lnmax(σrc,j′,pa,B)pa,B)−1,0], 29001_ 29010 24891 29034 12349 67482370 28996_ 29042 67273472 29009 8704 67273472 28955 560_ 29042 29027 24891 29034 24891 29040_ 29025 24891 28994 84054785 29040_ 29025 24891 28994 84054785 8704 28721 24891 28720 84267779 24891 (S39) where Dr 28996_ 29042 is expressed as a decimal. When the calculated confinement is below pa,B 29040_ 29025 24891 28994, this lower bound sets the logarithmic term to zero. The angle coefficient c 29027_ 28960 converts IR,j 29001_ 29010 24891 29034 into the local dilatancy angle: ψd,j=π180c IR,j. 28960_ 29028 24891 29034 12349 28953 28721 28728 28720 29027_ 28960 29001_ 29010 24891 29034 314 (S40) The factor π/180 28953 68408078 28721 28728 28720 converts the degree-valued intermediate c IR,j 29027_ 28960 29001_ 29010 24891 29034 to radians. A positive IR,j 29001_ 29010 24891 29034 therefore produces a positive dilation tendency, whereas IR,j=0 29001_ 29010 24891 29034 12349 28720 gives ψd,j=0 28960_ 29028 24891 29034 12349 28720. The dilation angle describes the tendency to expand, but the amount developed during one penetration step also depends on the interface displacement mobilized in that step. The current penetration increment is Δz=zn+1−zn 28673 29050 12349 29050_ 29038 8235 28721 8704 29050_ 29038. The parameter δy 28942_ 29049 denotes the characteristic mobilization displacement, and δcum,jn 28942_ 29027 29045 29037 24891 29034 29038 is the cumulative interfacial displacement committed through the preceding step. Equation (S41) combines these quantities to give the potential displacement for the current step: Δδpot,j=δy[1−exp(−ΔzUNKNOWNy)]exp(−δcum,jnUNKNOWNy). 28673 28942_ 29040 29039 29044 24891 29034 12349 28942_ 29049 67482370 28721 8704 67273472 8704 28673 29050 28942_ 29049 84054785 84267779 67273472 8704 28942_ 29027 29045 29037 24891 29034 29038 28942_ 29049 84054785 314 (S41) The bracketed term describes displacement mobilization over the current penetration increment. The final exponential factor reduces further mobilization as previously committed interface displacement accumulates. For Δz≥0 28673 29050 12821 28720, the resulting Δδpot,j 28673 28942_ 29040 29039 29044 24891 29034 is non-negative. S1.4.3 State Capacity and Void-Ratio Correction To determine the limit on further dilation, MD 29005_ 28996 compares the current void ratio with the critical-state value. This critical-state value is calculated from the mean effective stress in the converged MG 29005_ 28999–MS 29005_ 29011 state. When the current void ratio is lower, the positive difference is used to calculate the maximum dilation strain: pj′D,∗ 29040_ 29034 560 28996 24891 8707 =max[σ|,j′∗+2σh,j′∗3,pmin′], 12349 67482370 28955_ 69640972 24891 29034 560 8707 8235 28722 28955_ 29032 24891 29034 560 8707 28723 24891 29040 560_ 84267779 24891 (S42) ecs,jD,∗ 29029_ 29027 29043 24891 29034 28996 24891 8707 =eUNKNOWN−λc(pj′D,∗pa)UNKNOWNc, 12349 29029_ 28672 8704 28949_ 29027 67273472 29040_ 29034 560 28996 24891 8707 29040_ 29025 84054785 28948_ 29027 24891 ΨjD,∗ 28681_ 29034 28996 24891 8707 =ej∗−ecs,jD,∗, 12349 29029_ 29034 8707 8704 29029_ 29027 29043 24891 29034 28996 24891 8707 24891 Δε¯|,max,jD 28673 28962_ 69640972 24891 24891 29034 28996 =max(−ΨjD,∗,0)1+ej∗. 12349 67273472 8704 28681_ 29034 28996 24891 8707 24891 28720 84054785 28721 8235 29029_ 29034 8707 314 where ΨjD,∗<0 28681_ 29034 28996 24891 8707 12604 28720 indicates that the current void ratio is below the critical-state value. Dividing this positive void-ratio gap by 1+ej∗ 28721 8235 29029_ 29034 8707 gives the maximum cross-sectionally averaged dilation strain Δε¯|,max,jD 28673 28962_ 69640972 24891 24891 29034 28996. The current increment also depends on two distinct histories. The cumulative interface displacement δcum,jn 28942_ 29027 29045 29037 24891 29034 29038 in Eq. (S41) controls further displacement mobilization. The cumulative dilation-stage strain amplitude εcum,jD,n 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 records the dilation already applied through step n 29038. The factor ωj 28961_ 29034 reduces the current increment as this accumulated dilation approaches the current capacity: ωj=max[0,1−εcum,jD,nUNKNOWN"¯|,max,jD],UNKNOWN"¯|,max,jD>0andIR,j>0,UNKNOWN"¯|,max,jD≤0orIR,j=0. 28961_ 29034 12349 cases 67482370 28720 24891 28721 8704 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 28673 28962_ 69640972 24891 24891 29034 28996 84267779 24891& 28673 28962_ 69640972 24891 24891 29034 28996 12606 28720\ and\ 29001_ 29010 24891 29034 12606 28720 24891\\[12.0pt] 28720 24891& 28673 28962_ 69640972 24891 24891 29034 28996 12820 28720\ or\ 29001_ 29010 24891 29034 12349 28720 314 cases (S43) The local wall expansion is converted into an equivalent strain averaged over the caisson cross-section. The product Δδpot,jsinψd,j 28673 28942_ 29040 29039 29044 24891 29034 28960_ 29028 24891 29034 gives the outward displacement normal to the wall, and 4/Di 28724 68408078 28996_ 29033 is the perimeter-to-area ratio of the internal caisson section. Multiplying by ωj 28961_ 29034 gives the strain after history attenuation, before the remaining-capacity limit is applied: Δε¯|,raw,jD=4Δδpot,jsinψd,jDiωj. 28673 28962_ 69640972 24891 29042 29025 29047 24891 29034 28996 12349 28724 28673 28942_ 29040 29039 29044 24891 29034 28960_ 29028 24891 29034 28996_ 29033 28961_ 29034 314 (S44) The remaining capacity is obtained by subtracting the dilation already applied from the maximum amplitude: Δε¯|,rem,jD=max[Δε¯|,max,jD−εcum,jD,n,0]. 28673 28962_ 69640972 24891 29042 29029 29037 24891 29034 28996 12349 67482370 28673 28962_ 69640972 24891 24891 29034 28996 8704 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 24891 28720 84267779 314 (S45) The applied amplitude is the smaller of the untruncated value and this remaining capacity. When this amplitude is positive, the full potential interface displacement is recorded as Δδmob,j 28673 28942_ 29037 29039 29026 24891 29034; otherwise, the recorded displacement is zero: Δε¯|,jD 28673 28962_ 69640972 24891 29034 28996 =min[Δε¯|,raw,jD,Δε¯|,rem,jD], 12349 67482370 28673 28962_ 69640972 24891 29042 29025 29047 24891 29034 28996 24891 28673 28962_ 69640972 24891 29042 29029 29037 24891 29034 28996 84267779 24891 (S46) Δδmob,j 28673 28942_ 29037 29039 29026 24891 29034 =UNKNOWNUNKNOWNpot,j,UNKNOWN"¯|,jD>0,UNKNOWN"¯|,jD=0.. 12349 cases 28673 28942_ 29040 29039 29044 24891 29034 24891& 28673 28962_ 69640972 24891 29034 28996 12606 28720 24891\\[2.0pt] 28720 24891& 28673 28962_ 69640972 24891 29034 28996 12349 28720 314 cases 314 The applied dilation-stage strain amplitude gives the void-ratio increment: ΔejD=(1+ej∗)Δε¯|,jD. 28673 29029_ 29034 28996 12349 67273472 28721 8235 29029_ 29034 8707 84054785 28673 28962_ 69640972 24891 29034 28996 314 (S47) Adding this increment to the converged MS 29005_ 29011 void ratio and limiting the result to the admissible bounds gives the corrected void ratio for the current penetration step: ejn+1=clip[ej∗+ΔejD,emin,emax]. 29029_ 29034 29038 8235 28721 12349 29027 29036 29033 29040 67482370 29029_ 29034 8707 8235 28673 29029_ 29034 28996 24891 29029_ 24891 29029_ 84267779 314 (S48) The applied dilation-stage amplitude also updates the cumulative history carried into the subsequent penetration step: εcum,jD,n+1=εcum,jD,n+Δε¯|,jD. 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 8235 28721 12349 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 8235 28673 28962_ 69640972 24891 29034 28996 314 (S49) S1.5 Integrated Plug-Heave Assembly and Output Quantities The final void-ratio field returned by MD 29005_ 28996 is used to calculate plug length and heave by solid-volume conservation. The resulting geometry provides the domain for the refreshed hydraulic, stress, and pump-flow outputs. S1.5.1 Final Geometry, Output Refresh, and State Commitment The final geometry is assembled from the material layers reached by the caisson at zn+1 29050_ 29038 8235 28721. For material node j 29034, the active-node flag ℐjn+1 29001_ 29034 29038 8235 28721 equals 1 when ζ0,j≤zn+1 28944_ 28720 24891 29034 12820 29050_ 29038 8235 28721 and 0 otherwise. The active set collects the nodes with ℐjn+1=1 29001_ 29034 29038 8235 28721 12349 28721: n+1=j|ℐjn+1=1. 28993_ 29038 8235 28721 12349 \ 29034 12906 29001_ 29034 29038 8235 28721 12349 28721 \ 314 (S50) For each active material interval, solid-volume conservation converts its initial thickness Δζ0,j 28673 28944_ 28720 24891 29034 into its current contribution to plug length. The ratio (1+ejn+1)/(1+e0) 67273472 28721 8235 29029_ 29034 29038 8235 28721 84054785 68408078 67273472 28721 8235 29029_ 28720 84054785 accounts for the change in bulk soil volume at fixed solid volume. The area ratio αA 28939_ 28993 maps the reference layer volume associated with the outer caisson cross-section to the internal cross-section. Summing these contributions gives the plug length: Lplugn+1=∑j∈n+1αA1+ejn+11+e0Δζ0,j, 29004_ 29040 29036 29045 29031 29038 8235 28721 12349 4944 _ 29034 12850 28993_ 29038 8235 28721 28939_ 28993 28721 8235 29029_ 29034 29038 8235 28721 28721 8235 29029_ 28720 28673 28944_ 28720 24891 29034 24891 (S51) where Δζ0,j 28673 28944_ 28720 24891 29034 is the trapezoidal quadrature weight assigned to node j 29034 over the active material interval. The assembled plug length contains the penetrated portion below the seabed and the heaved portion above it. Subtracting the current penetration depth gives the net plug heave: hea|en+1=Lplugn+1−zn+1. 29032_ 29032 29029 29025 69640972 29029 29038 8235 28721 12349 29004_ 29040 29036 29045 29031 29038 8235 28721 8704 29050_ 29038 8235 28721 314 (S52) Equations (S51) and (S52) are the cellwise discretization of Eq. (S1). The corrected geometry also sets the domain height used to refresh the reported hydraulic and stress fields. The output-state height is reconstructed from the penetration depth and plug heave: Hn+1out=zn+1+hea|en+1. 29000_ 29038 8235 28721 29039 29045 29044 12349 29050_ 29038 8235 28721 8235 29032_ 29032 29029 29025 69640972 29029 29038 8235 28721 314 (S53) The resulting Hn+1out 29000_ 29038 8235 28721 29039 29045 29044 equals the assembled plug length. The superscript out 29039 29045 29044 marks its use in the output refresh. The final void-ratio field maps each active material coordinate to its output coordinate: xjout=∫0UNKNOWN0,jαA1+en+1(ζ0)1+e0dζ0, 29048_ 29034 29039 29045 29044 12349 4946 _ 28720 28944_ 28720 24891 29034 28939_ 28993 28721 8235 29029 29038 8235 28721 67273472 28944_ 28720 84054785 28721 8235 29029_ 28720\, 29028 28944_ 28720 24891 (S54) where en+1(ζ0) 29029 29038 8235 28721 67273472 28944_ 28720 84054785 denotes the final void-ratio field expressed over the material coordinate. The output seepage length is the finite-domain characteristic length evaluated with the refreshed plug length, ℓsout=ri[k|2krln(Re∗,outri)]1/2, 352_ 29043 29039 29045 29044 12349 29042_ 29033 67482370 29035_ 69640972 28722 29035_ 29042 67273472 29010_ 29029 8707 24891 29039 29045 29044 29042_ 29033 84054785 84267779 28721 68408078 28722 24891 (S55) where Re∗,out 29010_ 29029 8707 24891 29039 29045 29044 follows the adopted outer-domain rule in Section S1.3.1, evaluated with H=Hn+1out 29000 12349 29000_ 29038 8235 28721 29039 29045 29044. The refreshed hydraulic and vertical stress fields are evaluated from the same finite-domain relations: ijout=Δu(zn+1) wℓsoutsinh(Hn+1out/ℓsout)cosh(xjoutℓsout), 29033_ 29034 29039 29045 29044 12349 28673 29045 67273472 29050_ 29038 8235 28721 84054785 28941_ 29047 352_ 29043 29039 29045 29044 67273472 29000_ 29038 8235 28721 29039 29045 29044 68408078 352_ 29043 29039 29045 29044 84054785 67273472 29048_ 29034 29039 29045 29044 352_ 29043 29039 29045 29044 84054785 24891 (S56) σ|,j′out=γ′xjout−Δu(zn+1)sinh(xjout/ℓsout)sinh(Hn+1out/ℓsout). 28955_ 69640972 24891 29034 560 29039 29045 29044 12349 28941 560 29048_ 29034 29039 29045 29044 8704 28673 29045 67273472 29050_ 29038 8235 28721 84054785 67273472 29048_ 29034 29039 29045 29044 68408078 352_ 29043 29039 29045 29044 84054785 67273472 29000_ 29038 8235 28721 29039 29045 29044 68408078 352_ 29043 29039 29045 29044 84054785 314 (S57) The refreshed vertical stress is first bounded below by pmin′ 29040 560_ . The converged lateral stress is then projected onto the bounds associated with this refreshed vertical stress: σ^|,j′out=max[σ|,j′out,pmin′],σh,j′out=proj[σh,j′∗;pmin′,Kpσ^|,j′out]. 28955_ 69640972 24891 29034 560 29039 29045 29044 12349 67482370 28955_ 69640972 24891 29034 560 29039 29045 29044 24891 29040 560_ 84267779 24891 28955_ 29032 24891 29034 560 29039 29045 29044 12349 29040 29042 29039 29034 67482370 28955_ 29032 24891 29034 560 8707 24635 29040 560_ 24891 29003_ 29040 28955_ 69640972 24891 29034 560 29039 29045 29044 84267779 314 (S58) The refreshed mean and deviatoric stresses stored for output are pj′out 29040_ 29034 560 29039 29045 29044 =σ^|,j′out+2σh,j′out3, 12349 28955_ 69640972 24891 29034 560 29039 29045 29044 8235 28722 28955_ 29032 24891 29034 560 29039 29045 29044 28723 24891 (S59) p^j′out 29040_ 29034 560 29039 29045 29044 =max[pj′out,pmin′], 12349 67482370 29040_ 29034 560 29039 29045 29044 24891 29040 560_ 84267779 24891 qjout 29041_ 29034 29039 29045 29044 =|σ^|,j′out−σh,j′out|. 12349 69640972 28955_ 69640972 24891 29034 560 29039 29045 29044 8704 28955_ 29032 24891 29034 560 29039 29045 29044 69640972 314 After the MG 29005_ 28999–MS 29005_ 29011 iteration has ended, the solver evaluates this refresh once from the final corrected geometry and stores the resulting hydraulic and stress fields as outputs. At the end of the step, the final void ratio and applied interfacial displacement are committed: ej,n+1=ejn+1,δcum,jn+1=δcum,jn+Δδmob,j. 29029_ 29034 24891 29038 8235 28721 12349 29029_ 29034 29038 8235 28721 24891 28942_ 29027 29045 29037 24891 29034 29038 8235 28721 12349 28942_ 29027 29045 29037 24891 29034 29038 8235 28673 28942_ 29037 29039 29026 24891 29034 314 (S60) S1.5.2 Equivalent Pump-Flow Post-Processing After the final geometry and hydraulic field have been refreshed, the model calculates the equivalent pump flow associated with the prescribed suction. This quantity combines the Darcy inflow across the plug top with the rate of change in internal plug volume. The hydraulic quantities used below are taken from the refreshed output state: H=Hn+1out 29000 12349 29000_ 29038 8235 28721 29039 29045 29044, ℓs=ℓsout 352_ 29043 12349 352_ 29043 29039 29045 29044, and i=iout 29033 12349 29033 29039 29045 29044. The internal cross-sectional area is Ai=πDi24. 28993_ 29033 12349 28953 28996_ 29033 28722 28724 314 (S61) Using the refreshed hydraulic field, we calculate the upward flow rates across the plug top, Qin,top 29009_ 29033 29038 24891 29044 29039 29040, and the caisson tip, Qtip 29009_ 29044 29033 29040. Both are defined as positive upward: Qin,top=−Qx(0)=Aik|i(0,z)=Aik|Δu(z) wℓsinh(H/ℓs),Qtip=−Qx(H)=Aik|itip. 29009_ 29033 29038 24891 29044 29039 29040 12349 8704 29009_ 29048 67273472 28720 84054785 12349 28993_ 29033 29035_ 69640972 29033 67273472 28720 24891 29050 84054785 12349 28993_ 29033 29035_ 69640972 28673 29045 67273472 29050 84054785 28941_ 29047 352_ 29043 (H/ _s) 24891 29009_ 29044 29033 29040 12349 8704 29009_ 29048 67273472 29000 84054785 12349 28993_ 29033 29035_ 69640972 29033_ 29044 29033 29040 314 (S62) The second contribution is the change in internal plug volume during the penetration step. Because H=z+hea|e 29000 12349 29050 8235 29032_ 29032 29029 29025 69640972 29029, this change contains both the penetration increment and the heave increment. For successive accepted steps n 29038 and n+1 29038 8235 28721, let Δz=zn+1−zn 28673 29050 12349 29050_ 29038 8235 28721 8704 29050_ 29038, let z˙ 29050 denote the penetration rate, and let Δt 28673 29044 denote the elapsed time: Qgeom=Ai(zn+1−znUNKNOWNt+hea|en+1−hea|enUNKNOWNt),Δt=Δz˙. 29009_ 29031 29029 29039 29037 12349 28993_ 29033 67273472 29050_ 29038 8235 28721 8704 29050_ 29038 28673 29044 8235 29032_ 29032 29029 29025 69640972 29029 29038 8235 28721 8704 29032_ 29032 29029 29025 69640972 29029 29038 28673 29044 84054785 24891 28673 29044 12349 28673 29050 29050 314 (S63) The equivalent pump-flow estimate is the sum of the top inflow and the geometric volume-change rate: Qpump=Qin,top+Qgeom. 29009_ 29040 29045 29037 29040 12349 29009_ 29033 29038 24891 29044 29039 29040 8235 29009_ 29031 29029 29039 29037 314 (S64) The solver records Qpump 29009_ 29040 29045 29037 29040 after the output refresh. Calculating the physical flow rate requires a measured time interval or penetration rate. For benchmark profiles without timing data, the calculation uses a unit time interval, Δt=1 28673 29044 12349 28721; the resulting pump-flow values are excluded from the plug-heave benchmark assessment. S1.6 Discrete Solution and State Transfer The numerical solver advances through the prescribed penetration depths, using the accepted state from one step to initialize the next (Supplementary Fig. 1). Supplementary Fig. 1. Stepwise numerical solution and state transfer of the soil-plug model. Trial state denotes fields updated during the MG 29005_ 28999–MS 29005_ 29011 iteration; converged state denotes the accepted solution before the MD 29005_ 28996 correction. Committed state carries the accepted state and history into the next penetration step. Output state contains the post-MD 29005_ 28996 refreshed fields stored for reporting and verification. S1.6.1 Material Grid and Pre-Step State The solver represents the initial soil profile on a fixed material grid. The coordinate ζ0,j 28944_ 28720 24891 29034 is the initial depth of node j 29034 below the seabed, with j=1,…,N 29034 12349 28721 24891… 24891 29006. At each penetration step, the solver applies the active-node rule to this grid. Active nodes contribute to the plug-heave integral and receive the MS 29005_ 29011 and MD 29005_ 28996 updates. For nodes that were already active at the accepted step n 29038, the solver reads the following committed nodal state before calculating step n+1 29038 8235 28721: ej,n,σ|,j,n′H,σh,j,n′H,pj,n′H,η0,j,μ¯j,n,δcum,jn,εcum,jD,n. \ 29029_ 29034 24891 29038 24891\, 28955_ 69640972 24891 29034 24891 29038 560 29000 24891\, 28955_ 29032 24891 29034 24891 29038 560 29000 24891\, 29040_ 29034 24891 29038 560 29000 24891\, 28945_ 28720 24891 29034 24891\, 28950_ 29034 24891 29038 24891\, 28942_ 29027 29045 29037 24891 29034 29038 24891\, 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 \ 314 (S65) The superscript H 29000 marks stress quantities carried from the preceding accepted step. The void ratio ej,n 29029_ 29034 24891 29038 is the final value after the preceding MD 29005_ 28996 correction. The vertical stress σ|,j,n′H 28955_ 69640972 24891 29034 24891 29038 560 29000 is obtained from the refreshed final geometry, while σh,j,n′H 28955_ 29032 24891 29034 24891 29038 560 29000 is the lateral stress synchronized at MG 29005_ 28999–MS 29005_ 29011 convergence. The rebound calculation uses the positive reference stress pj,n′H=max[cp(K0)σ|,j,n′H,pmin′] 29040_ 29034 24891 29038 560 29000 12349 67482370 29027_ 29040 67273472 29003_ 28720 84054785 28955_ 69640972 24891 29034 24891 29038 560 29000 24891 29040 560_ 84267779, where cp(K0)=(1+2K0)/3 29027_ 29040 67273472 29003_ 28720 84054785 12349 67273472 28721 8235 28722 29003_ 28720 84054785 68408078 28723. The remaining quantities record the constitutive history of each node. The reference η0,j 28945_ 28720 24891 29034 stores the stress ratio at first activation, and μ¯j,n 28950_ 29034 24891 29038 stores the largest path mobilization reached through step n 29038. The variables δcum,jn 28942_ 29027 29045 29037 24891 29034 29038 and εcum,jD,n 28962_ 29027 29045 29037 24891 29034 28996 24891 29038 record the applied MD 29005_ 28996 history. The trial calculations use this committed state as a fixed reference. After the current step is accepted, its updated values become the pre-step state for the next depth. A node that first becomes active at step n+1 29038 8235 28721 starts from the initial void ratio e0 29029_ 28720, its geostatic vertical stress γ′ζ0,j 28941 560 28944_ 28720 24891 29034, and the corresponding K0 29003_ 28720 lateral stress. Its accumulated path mobilization, interfacial displacement, and dilation-stage strain amplitude start from zero. The activation reference η0,j 28945_ 28720 24891 29034 is recorded from the synchronized state obtained during this first active step. In the stress relations, K0 29003_ 28720 denotes the at-rest earth-pressure coefficient; K0(λ1ri) 29003_ 28720 67273472 28949_ 28721 29042_ 29033 84054785 denotes the modified Bessel function used in the hydraulic formulation. Each benchmark profile is represented by 121 uniformly spaced material nodes extending from the initial seabed to the final penetration depth. The input depths are sorted and duplicate depths are removed to define the penetration steps; each retained depth uses its paired suction value. Spatial volume integrals use trapezoidal quadrature. When a penetration depth falls between adjacent grid nodes, the terminal value is linearly interpolated and the remaining partial interval is included in the integral. S1.6.2 Fixed-Point Iteration and Convergence Each iteration begins with a plug-heave estimate hest(k) 29032_ 29029 29043 29044 67273472 29035 84054785 and a trial void-ratio field etrial(k) 29029_ 29044 29042 29033 29025 29036 67273472 29035 84054785. These values determine the plug geometry, seepage length, hydraulic field, and effective stresses. The stress update produces a candidate void-ratio field ecand(k) 29029_ 29027 29025 29038 29028 67273472 29035 84054785, which is integrated to obtain the candidate height hcand(k) 29032_ 29027 29025 29038 29028 67273472 29035 84054785. With the committed history from step n 29038 held fixed, this calculation defines the fixed-point map: (ecand(k),hcand(k))=ℱn+1(etrial(k),hest(k)),(e∗,hea|e∗)=ℱn+1(e∗,hea|e∗). 67273472 29029_ 29027 29025 29038 29028 67273472 29035 84054785 24891 29032_ 29027 29025 29038 29028 67273472 29035 84054785 84054785 12349 28998_ 29038 8235 28721\! 67273472 29029_ 29044 29042 29033 29025 29036 67273472 29035 84054785 24891 29032_ 29029 29043 29044 67273472 29035 84054785 84054785 24891 67273472 29029 8707 24891 29032_ 29032 29029 29025 69640972 29029 8707 84054785 12349 28998_ 29038 8235 28721\! 67273472 29029 8707 24891 29032_ 29032 29029 29025 69640972 29029 8707 84054785 314 (S66) At a fixed point, the input estimate and returned candidate describe the same coupled geometry and void-ratio state. The solver measures this closure through the difference between the height used to construct the trial fields and the height returned by the updated void ratio: rh(k)=|hcand(k)−hest(k)|. 29042_ 29032 67273472 29035 84054785 12349 69640972 29032_ 29027 29025 29038 29028 67273472 29035 84054785 8704 29032_ 29029 29043 29044 67273472 29035 84054785 69640972 314 (S67) The iteration is accepted when either the relative residual is below 10−5 28721 28720 8704 28725 or the absolute residual is below 10−8m 28721 28720 8704 28728\, 29037: rh(k)max(|hest(k)|,10−10m)<10−5orrh(k)<10−8m. 29042_ 29032 67273472 29035 84054785 67273472 69640972 29032_ 29029 29043 29044 67273472 29035 84054785 69640972 24891 28721 28720 8704 28721 28720\, 29037 84054785 12604 28721 28720 8704 28725 29042_ 29032 67273472 29035 84054785 12604 28721 28720 8704 28728\, 29037 314 (S68) The relative criterion scales the change by the current heave estimate. The absolute criterion provides a direct tolerance when the heave is close to zero, and the 10−10m 28721 28720 8704 28721 28720\, 29037 denominator bound prevents division by zero. When the criterion is not satisfied, the next estimate is hest(k+1)=hest(k)+0.7[hcand(k)−hest(k)] 29032_ 29029 29043 29044 67273472 29035 8235 28721 84054785 12349 29032_ 29029 29043 29044 67273472 29035 84054785 8235 28720 314 28727 67482370 29032_ 29027 29025 29038 29028 67273472 29035 84054785 8704 29032_ 29029 29043 29044 67273472 29035 84054785 84267779, and the new void-ratio candidate becomes the trial field for the next iteration. This partial update limits oscillation in the feedback between geometry and soil state. Each penetration step is allowed at most 120 iterations. If no iteration satisfies Eq. (S68), the solver reports non-convergence and ends the step before state transfer. After the criterion is satisfied, the solver performs one synchronization pass at the accepted candidate geometry. This pass recalculates the MG 29005_ 28999–MS 29005_ 29011 fields so that the returned void ratio, stresses, and plug heave refer to the same geometry. The accepted iteration is denoted by kcon| 29035_ 29027 29039 29038 69640972. Its mobilization degree μj,n+1con| 28950_ 29034 24891 29038 8235 28721 29027 29039 29038 69640972 updates the historical maximum at node j 29034: μ¯j,n+1=max[μ¯j,n,μj,n+1con|]. 28950_ 29034 24891 29038 8235 28721 12349 67482370 28950_ 29034 24891 29038 24891 28950_ 29034 24891 29038 8235 28721 29027 29039 29038 69640972 84267779 314 (S69) The synchronized void-ratio field becomes the pre-dilation value passed to MD 29005_ 28996: ej,n+1∗=ej,n+1∗,(kcon|). 29029_ 29034 24891 29038 8235 28721 8707 12349 29029_ 29034 24891 29038 8235 28721 8707 24891 67273472 29035_ 29027 29039 29038 69640972 84054785 314 (S70) The same synchronization pass returns the corresponding vertical stress, lateral stress, mean effective stress, critical-state void ratio, and state parameter: σ|,j,n+1′∗ 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 8707 =σ|,j,n+1′(kcon|),σh,j,n+1′∗=σh,j,n+1′(kcon|), 12349 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 67273472 29035_ 29027 29039 29038 69640972 84054785 24891 28955_ 29032 24891 29034 24891 29038 8235 28721 560 8707 12349 28955_ 29032 24891 29034 24891 29038 8235 28721 560 67273472 29035_ 29027 29039 29038 69640972 84054785 24891 (S71) pj,n+1′∗ 29040_ 29034 24891 29038 8235 28721 560 8707 =pj,n+1′(kcon|),ecs,j,n+1∗ 12349 29040_ 29034 24891 29038 8235 28721 560 67273472 29035_ 29027 29039 29038 69640972 84054785 24891 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 8707 =ecs,j,n+1(kcon|),Ψj,n+1∗ 12349 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 67273472 29035_ 29027 29039 29038 69640972 84054785 24891 28681_ 29034 24891 29038 8235 28721 8707 =Ψj,n+1∗,(kcon|). 12349 28681_ 29034 24891 29038 8235 28721 8707 24891 67273472 29035_ 29027 29039 29038 69640972 84054785 314 These quantities form the synchronized pre-dilation state: j,n+1∗=ej,n+1∗,σ|,j,n+1′∗,σh,j,n+1′∗,pj,n+1′∗,ecs,j,n+1∗,Ψj,n+1∗,μ¯j,n+1. 29011_ 29034 24891 29038 8235 28721 8707 12349 \ 29029_ 29034 24891 29038 8235 28721 8707 24891 28955_ 69640972 24891 29034 24891 29038 8235 28721 560 8707 24891 28955_ 29032 24891 29034 24891 29038 8235 28721 560 8707 24891 29040_ 29034 24891 29038 8235 28721 560 8707 24891 29029_ 29027 29043 24891 29034 24891 29038 8235 28721 8707 24891 28681_ 29034 24891 29038 8235 28721 8707 24891 28950_ 29034 24891 29038 8235 28721 \ 314 (S72) The star marks quantities evaluated in this synchronized pre-dilation state. MD 29005_ 28996 acts once on this state, after which the corrected void ratio determines the final geometry and output refresh. The accepted histories are then committed for the next penetration depth. S1.7 Contribution and Update-Order Sensitivity of MD 29005_ 28996 To evaluate the influence of MD 29005_ 28996 after MG 29005_ 28999–MS 29005_ 29011 convergence, we performed two numerical checks across the 14 benchmark profiles. One compared predictions from MG 29005_ 28999 alone, MG 29005_ 28999–MS 29005_ 29011, and the complete MG 29005_ 28999–MS 29005_ 29011–MD 29005_ 28996 calculation to quantify the contribution of MD 29005_ 28996. The other tested whether one MD 29005_ 28996 update per penetration step was sufficient by comparing it with repeated updates at the same penetration depth. In the reported calculation, each penetration step completed MG 29005_ 28999–MS 29005_ 29011 convergence before MD 29005_ 28996 corrected the void-ratio field and plug heave. After this single correction, the calculation advanced to the next penetration depth. Thus, “once” means one MD 29005_ 28996 correction per penetration step. Because this correction also changes the plug geometry used by MG 29005_ 28999–MS 29005_ 29011, the comparison calculation held the penetration depth fixed and repeated the update. The MD 29005_ 28996-corrected heave and void ratio were returned to MG 29005_ 28999–MS 29005_ 29011, after which MD 29005_ 28996 was applied again. Repetition stopped when two successive post-MD 29005_ 28996 heaves differed by less than 10−5 28721 28720 8704 28725 m, with a maximum of eight repetitions. The nine final-state profiles comprised four Kim and Kim cases (T1-1 S 25g, T1-2 S 50g, T3-1 S 50g, and T4-3 S 25g) and five Tran and Randolph cases (TRAN-T1–TRAN-T5). The five continuous-process profiles comprised Tjelta TJELTA-A3, Wang full process, and Zhao ZHAO-T1–ZHAO-T3. Supplementary Table 3. Contribution and update-order sensitivity of MD 29005_ 28996 across the 14 benchmark profiles. Check Metric Value Component combination Mean absolute final-heave error with MG 29005_ 28999 only 58.42% Component combination Mean absolute final-heave error with MG 29005_ 28999 and MS 29005_ 29011 9.01% Component combination Mean absolute final-heave error with MG 29005_ 28999, MS 29005_ 29011, and MD 29005_ 28996 9.48% Relative MD 29005_ 28996 heave increment Mean absolute ratio of the MD 29005_ 28996 heave increment to final plug heave, all profiles 0.956% Relative MD 29005_ 28996 heave increment Mean absolute ratio of the MD 29005_ 28996 heave increment to final plug heave, final state prediction 0.097% Relative MD 29005_ 28996 heave increment Mean absolute ratio of the MD 29005_ 28996 heave increment to final plug heave, continuous process prediction 2.503% Relative MD 29005_ 28996 heave increment Maximum absolute ratio of the MD 29005_ 28996 heave increment to final plug heave 3.478% (Zhao T3) Absolute MD 29005_ 28996 heave increment Mean absolute MD 29005_ 28996 heave increment, all profiles 0.0675 cm Absolute MD 29005_ 28996 heave increment Mean absolute MD 29005_ 28996 heave increment, final state prediction 0.0475 cm Absolute MD 29005_ 28996 heave increment Mean absolute MD 29005_ 28996 heave increment, continuous process prediction 0.1035 cm Single versus repeated MD 29005_ 28996 update Mean absolute final-heave difference, all profiles 0.0091 cm Single versus repeated MD 29005_ 28996 update Maximum absolute final-heave difference 0.0200 cm (Wang full process) Single versus repeated MD 29005_ 28996 update Maximum absolute difference between final-heave errors 0.418% (Wang full process) Adding MS 29005_ 29011 to MG 29005_ 28999 reduced the mean absolute final-heave error from 58.42% to 9.01%, representing the dominant change in the component comparison. Including MD 29005_ 28996 changed the aggregate error to 9.48%, while its mean absolute contribution was 0.0675 cm, or 0.956% of the final heave. This contribution was more visible in the continuous-process profiles than in the final-state profiles, averaging 2.503% and 0.097%, respectively, and reaching 3.478% for Zhao T3. These results identify MD 29005_ 28996 as a modest, case-dependent correction to the predicted heave; the principal reduction in final-heave error is associated with MS 29005_ 29011. Repeating the MD 29005_ 28996 update at the same penetration depth changed the final heave by 0.0091 cm on average and by at most 0.0200 cm. The largest absolute difference between the two final-heave errors was 0.418%, also for the Wang full-process profile. The close agreement across the 14 profiles supports one MD 29005_ 28996 correction after MG 29005_ 28999–MS 29005_ 29011 convergence in each penetration step. S1.8 Numerical Safeguards and Model Scope Three safeguards keep the calculation well defined near limiting states. The solver uses pmin′=0.1 29040 560_ 12349 28720 314 28721 kPa as a positive stress floor and ϵUNKNOWN=10−8 28943_ 28945 12349 28721 28720 8704 28728 as the lower bound in the mobilization denominator. The void-ratio update is limited to the case-specific, source-based bounds emin 29029_ and emax 29029_ , which define the admissible density range. The bounded stresses enter the constitutive update, but the solver also records the unbounded vertical stress. A value of σ|′≤0 28955 560_ 69640972 12820 28720 indicates that the assumed seepage field has reduced the computed vertical effective stress to zero or below. The hydraulic equations are evaluated only when H>0 29000 12606 28720. At the initial state H=0 29000 12349 28720, the zero-suction branch assigns zero hydraulic gradients and flow rates. The seepage model assumes saturated sand and Darcy flow, with quasi-steady conditions over each penetration step. It uses the equivalent outer-domain influence radius Re∗ 29010_ 29029 8707 and holds k| 29035_ 69640972 and kr 29035_ 29042 constant within that step. The model does not resolve piping or particle migration. To monitor this boundary, the solver retains the hydraulic gradient at the caisson tip, itip 29033_ 29044 29033 29040, at every penetration step. The recorded itip 29033_ 29044 29033 29040 can be compared with a piping criterion appropriate to the soil and hydraulic boundary. A step that reaches or exceeds that criterion lies at the limit of the present continuum calculation. The formulation is intended for saturated sand when seepage can be represented by Darcy flow and the soil–caisson interface response can be represented by a cross-sectional average. The model would require extension where permeability changes strongly during installation, deformation localizes into shear bands, or the interface response becomes strongly three-dimensional. Complex elastoplastic stress paths and soil-fabric reconstruction after failure or liquefaction also lie outside its present scope. Supplementary Note 2. Workflow Audit Procedures and Results The main text reports the results of two audits that examined the equations, code, and check records retained during model development. This note explains what each audit role received, which answers were withheld, how its response was scored, and what problems it found. S2.1 Materials Shown to Each Audit Role We used two audit tasks to test whether an audit role could find problems without being given their answers. In the blinded replay, the role received materials related to a problem previously identified during development, but not its label or answer key. The audit role then had to identify and locate the problem from those materials. For each derivation-to-code probe, a code-writing role first implemented the supplied equations and completed a fixed set of checks. A separate audit role then examined the equations, generated code, and check results to identify problems that those checks did not cover. Supplementary Table 4 specifies the materials shown and withheld in each task. Supplementary Table 4. Materials shown to and withheld from each audit role. Audit task Materials shown to the audit role Information not shown Question evaluated Blinded audit replay (A1–A9) The relevant physical relation, calculation rule, required input and output of the code module, uncorrected code or output, and available check records The issue label, answer key, prior discussion, correction history, and other project files Can the audit role find and locate a known issue without being given its answer? Derivation-to-code probes (P1–P5) The supplied equations and state-update rules, generated code, and results of the fixed implementation checks The target problem and answer Can the audit role identify omissions or inconsistencies not covered by the predeclared checks? S2.2 How the Audits Were Scored We scored both audit tasks using criteria that linked each reported problem to a specific equation, calculation step, or code location. For the blinded replay, we compared each report with an answer key held outside the audit packet. A replay passed when it identified the target problem and its location, explained the physical or numerical consequence, and proposed a valid correction with a follow-up check. We recorded false positives and human assistance separately. For each derivation-to-code probe, a report passed when it identified a specific disagreement between the supplied equations or update rules and the generated code, located the affected code, and cited the supporting check result. Supplementary Table 5 gives the complete scoring criteria. Supplementary Table 5. What counted as a pass in each audit. Scored item What was checked Criterion for a pass Finding and location Whether the report found the target problem and named the affected equation, calculation step, or code location The target problem and its location are both identified correctly Consequence and correction Physical or numerical consequence, proposed correction, and follow-up check The correction addresses the problem and the follow-up check can verify it False positives and human assistance Any unrelated finding and any help given during replay The target problem is found without an incorrect additional finding or an extra hint Independent probe finding The equation or rule that the code failed to follow, its code location, and the supporting check result The report identifies a specific omission or disagreement from the supplied equations, code, and check results S2.3 Blinded Audit-Replay Results The audit role found and located all nine target problems without false positives or human assistance. Supplementary Table 6 lists the target problem in each replay and the check needed to identify it. Supplementary Table 6. Blinded audit-replay results. Case Target problem What the audit role had to check A1 Correction for near-wall soil dilation absent from the solver call chain Whether the solver calls the MD 29005_ 28996 correction A2 Degree-valued ψd 28960_ 29028 passed directly to a trigonometric function expecting radians Whether the angle unit matches the function input A3 Committed local increment not limited by the remaining dilation capacity Whether the applied increment remains within the available capacity A4 The reported mean effective stress and deviatoric stress retained their trial values and were not recalculated from the final stresses Whether the reported mean effective stress and deviatoric stress are recalculated from the final stresses A5 Trial iterations overwrite committed history before convergence is accepted Whether history is saved only after convergence A6 Reference confinement uses vertical trial stress instead of lateral trial stress Which stress is used as the reference confinement A7 Pump-flow post-processing uses the tip-flow component instead of the top inflow component Which flow component enters the reported pump flow A8 Final hydraulic state reuses the trial values instead of recomputing from final geometry Whether the hydraulic field is recalculated from the final geometry A9 Post-rebound capacity is computed from the pre-rebound void ratio instead of the rebound-updated value Which void ratio is used to calculate the remaining capacity For A6, the initial audit packet omitted the stress-regularization relation needed to determine the reference confinement. We supplied this equation before scoring. The issue label and answer key remained withheld. S2.4 Derivation-to-Code Probe Results All 36 fixed implementation checks passed across the five probes, yet the independent audit identified five additional problems in P1–P4 (Supplementary Table 7). The findings involved code-branch handling, incomplete instructions, an incorrect count of checked trials, and a disagreement between an equation and its implementation. The audit found no additional problem when it repeated the final hydraulic and stress checks in P5. Supplementary Table 7. Derivation-to-code probe results. Probe Calculation checked Fixed checks passed What the independent audit found P1 Recalculation of final outputs in converged, edge, and nonconverged branches 10/10 One branch omitted the final-output recalculation; an additional check was added P2 Active nodes and the piecewise ω 28961 calculation in MD 29005_ 28996 11/11 One code branch could fail, and the supplied instructions did not state clearly whether inactive grid nodes remained in the input P3 Saving MS 29005_ 29011 history only after convergence is accepted 5/5 The code stopped checking trials as soon as it found one that met the convergence criterion. However, the recorded count included all supplied trials, including those that had not been checked. P4 Applying MS 29005_ 29011 rebound before calculating the remaining capacity 5/5 The code and equation treated Δereb 28673 29029_ 29042 29029 29026 differently in the numerator P5 Recalculating hydraulic and stress outputs from the final geometry 5/5 All final-refresh checks were repeated and no additional problem was found S2.5 What the Workflow Audits Establish Within the evaluated tasks, the replay and probe results support two distinct audit functions: recovering known problems from retained records and identifying implementation problems missed by fixed checks. Benchmark comparisons provide separate evidence for prediction accuracy. Supplementary Note 3. Checks of the Final Code After addressing the audit findings, we checked whether the final penetration-step solver followed the equations, update order, and state-commitment rules defined in Supplementary Note 1. Module-level tests examined individual calculations and state updates. Analytical reference checks compared selected relations with independent evaluations, whereas regression tests checked whether unchanged input sequences reproduced the stored outputs. S3.1 Module-Level Tests All 33 module-level tests passed. Supplementary Table 8 groups the tests by the part of the calculation examined and reports the main operations, acceptance criteria, and outcome for each group. Supplementary Table 8. Module-level tests and results. Test group Number of tests Main operations checked Acceptance criterion and observed result Outcome Hydraulic closure and moving-boundary stress 10 Pressure-difference recovery, tip gradient, free-surface stress, suction-controlled boundaries, and pump-flow separation Integral and boundary residuals remain within the adopted tolerance; pump flow does not update plug heave or void ratio Passed Bolton relationship and MD 29005_ 28996 closure 9 Node activation, reference confinement, dilatancy, mobilized displacement, capacity truncation, and final clipping The correction is non-negative, capacity-limited, and confined to active nodes Passed MS 29005_ 29011 path update and projection constraints 6 Mean effective stress, lateral projection, critical-state variables, rebound, state-gap release, and clipping Updated variables match analytical or independently evaluated values Passed State commitment and final output refresh 8 Separation of trial and accepted states, saving after convergence, final MD 29005_ 28996 update, output recalculation, and saved iteration and check results Histories are saved only from accepted states, and final outputs are recalculated from the final geometry and stresses Passed S3.2 Analytical Reference and Regression Checks All five analytical reference checks and four regression tests passed (Supplementary Table 9). Supplementary Table 9. Analytical reference and regression checks. ID Check Quantity or comparison Criterion and observed result Outcome R-01 Pressure-difference closure γw∫0Hi(x,z)x−Δu(z) 28941_ 29047 4946 _ 28720 29000 29033 67273472 29048 24891 29050 84054785\, 29028 29048 8704 28673 29045 67273472 29050 84054785 Absolute residual remains within the adopted tolerance Passed R-02 Moving free-surface stress σ|′(0,z) 28955 560_ 69640972 67273472 28720 24891 29050 84054785 Zero within numerical tolerance Passed R-03 Deep-limit approximation Finite-domain seepage kernel and its exponential deep-limit form Relative difference remains within tolerance Passed R-04 Critical-state consistency ecs 29029_ 29027 29043 and Ψ=e−ecs 28681 12349 29029 8704 29029_ 29027 29043 Matches independent scalar evaluation Passed R-05 MD 29005_ 28996 capacity truncation Δε¯|D≤Δε¯|,remD 28673 28962_ 69640972 28996 12820 28673 28962_ 69640972 24891 29042 29029 29037 28996 No node exceeds the remaining capacity Passed B-01 Input-sequence consistency Penetration depths, suction-pressure history, geometry, and soil inputs Input arrays and parameter assignments remain unchanged Passed B-02 MG 29005_ 28999–MS 29005_ 29011 baseline hea|eGS 29032_ 29032 29029 29025 69640972 29029 28999 29011, e∗ 29029 8707, and convergence status Stored main-chain outputs are reproduced Passed B-03 MD 29005_ 28996 contribution baseline Δhea|eD 28673 29032_ 29032 29029 29025 69640972 29029 28996, ΔeD 28673 29029 28996, and capacity flags Stored correction outputs are reproduced Passed B-04 Post-MD 29005_ 28996 refresh baseline Final plug heave, refreshed stresses, and pump-flow outputs Stored final outputs are reproduced Passed S3.3 Results of the Final-Code Checks Within the tested scope, these results show that the code follows the declared equations, update order, and state-commitment rules. Supplementary Note 4. Benchmark Data, Input Sources, and Comparison Results S4.1 Public Benchmark Files The public files separate the values used in the benchmark comparisons from the sources and uses of the model inputs. Supplementary Data 1 provides the case-level values for three benchmark comparisons: formula reduction and volume closure, final state prediction, and continuous process prediction. Supplementary Data 2 records each adopted global and case-specific input, its source and source class, and where it enters the calculation. Definitions of all fields in the two data files are provided in SUPPLEMENTARY_DATA_README.md. Supplementary Data 2 uses five source classes for adopted inputs: 1. direct: a value reported by a source. 2. derived: a value calculated from reported values. 3. design_parameter: a value adopted from site or design information. 4. analog: a value transferred from a cited related source. 5. assumption: a value fixed by a declared modeling rule. The label unavailable marks a quantity for which no sourced value was adopted. The label mixed marks a ledger row that combines adopted values from different source classes; the source-basis field identifies the class of each value. S4.2 Shared Input Rules The final-state and continuous-process predictions use common assignment rules for three model inputs. For eUNKNOWN/emax 29029_ 28672 68408078 29029_ , the central value is 1.0; varying this ratio from 0.9 to 1.1 generates the prediction sensitivity band. Published sand critical-state and state-parameter data support this range [48, 49, 50, 51]. Poisson’s ratio is fixed at ν=0.30 28951 12349 28720 314 28723 28720, within reported ranges for granular soils[37, 47]. Kim and Kim (2020)[16] report case-specific ki/kout 29035_ 29033 68408078 29035_ 29039 29045 29044 values ranging from 1.26 to 1.46; a representative value of 1.3 is adopted for those cases. The other prediction cases use ki/kout=3.0 29035_ 29033 68408078 29035_ 29039 29045 29044 12349 28723 314 28720, following suction-caisson installation guidance[6]. The equivalent model ratio kr/k| 29035_ 29042 68408078 29035_ 69640972 is assigned these same values as a model assumption. The formula reduction and volume closure cases are evaluated independently of these three inputs. S4.3 Outer-Domain Boundary Rule for Re∗ 29010_ 29029 8707 The outer-domain boundary rule determines whether the seepage kernel uses the unbounded modal estimate Re,∞ 29010_ 29029 24891 561 or a finite radius derived from the experimental boundary. Section S1.3.1 defines the adopted radius Re∗ 29010_ 29029 8707; when a source-based radius Rb 29010_ 29026 is available, it caps Re,∞ 29010_ 29029 24891 561. Because Rb 29010_ 29026 represents the experimental boundary, it remains fixed during the sensitivity analysis of eUNKNOWN/emax 29029_ 28672 68408078 29029_ . Supplementary Table 10. Source-based outer-domain boundary rules. Case group External-boundary basis Adopted rule Source class Kim and Kim (2020)[16] 900-m-diameter centrifuge container at 25g and 50g Rb=0.45Ng 29010_ 29026 12349 28720 314 28724 28725 29006_ 29031 m derived Tran and Randolph (2008)[14] 390-m by 650-m strongbox at 100g Rb=0.195×100=19.5 29010_ 29026 12349 28720 314 28721 28729 28725 8706 28721 28720 28720 12349 28721 28729 314 28725 m derived Tjelta (1995)[41] Full-scale field bucket without a test-box boundary Re∗=Re,∞ 29010_ 29029 8707 12349 29010_ 29029 24891 561; no finite Rb 29010_ 29026 direct Wang et al. (2019)[44] 1.50-m by 0.96-m by 0.90-m model box Rb=0.48 29010_ 29026 12349 28720 314 28724 28728 m derived Zhao et al. (2025)[43] 1.0-m by 0.5-m by 0.6-m model tank Rb=0.25 29010_ 29026 12349 28720 314 28722 28725 m derived Here Ng 29006_ 29031 is the centrifuge acceleration factor. For each prediction case with a finite cutoff, Re,∞ 29010_ 29029 24891 561 remains below that cutoff; Tjelta has no finite cutoff. The predictions therefore use Re,∞ 29010_ 29029 24891 561 in all evaluated prediction cases. S4.4 Parameter Sources The parameter-source ledger contains 122 entries. Supplementary Table 11 summarizes model inputs shared across runs, case-assignment rules, and numerical settings. Supplementary Table 12 shows how data obtained from individual studies enter each benchmark comparison. Supplementary Table 11. Shared model inputs, assignment rules, and numerical settings. Setting Adopted value or rule Source basis Use eUNKNOWN/emax 29029_ 28672 68408078 29029_ 1.0 central; 0.9–1.1 band Declared sensitivity rule Prediction sensitivity band λc 28949_ 29027 0.019 Same-form silica-sand CSL prior State update κc 28948_ 29027 0.7 Same-form silica-sand CSL prior State update pa 29040_ 29025 100 kPa Reference pressure paired with the CSL prior State update κs 28948_ 29043 0.006 Fixed global coefficient MS 29005_ 29011 update ν 28951 0.30 Granular-soil literature range Lateral-stress update kr/k| 29035_ 29042 68408078 29035_ 69640972 1.3 for Kim and Kim; 3.0 otherwise Model assumption: same value as the adopted ki/kout 29035_ 29033 68408078 29035_ 29039 29045 29044 Seepage kernel pmin′ 29040 560_ 0.1 kPa Numerical lower bound Stress and state update ϵUNKNOWN 28943_ 28945 10−8 28721 28720 8704 28728 Denominator lower bound MS 29005_ 29011 mobilization Q 29009 10 Bolton-type quartz/silica-sand prior MD 29005_ 28996 correction cUNKNOWN 29027_ 28960 0.5 Fixed Bolton-type coefficient MD 29005_ 28996 correction pa,B 29040_ 29025 24891 28994 1 kPa Reference-pressure convention MD 29005_ 28996 correction δy 28942_ 29049 0.001 m Fixed characteristic displacement Eq. (S41) Material grid 121 uniformly spaced nodes Fixed numerical configuration All solver runs Penetration sequence Sorted unique benchmark depths Direct benchmark input All solver runs Suction sequence Value paired with each retained depth; no solver interpolation Direct benchmark input All solver runs Spatial integration Trapezoidal rule with linear terminal-interval interpolation Fixed numerical configuration Plug-volume assembly Relaxation factor 0.7 Fixed numerical configuration Heave iteration Convergence Relative tolerance 10−5 28721 28720 8704 28725 or residual 10−8 28721 28720 8704 28728 m Eq. (S68) Heave iteration Iteration limit 120 per step; reject non-converged steps Fixed numerical configuration Run acceptance Supplementary Table 12. Study data and their use in each benchmark comparison. Study Benchmark comparison Data obtained from the study Use in the comparison Koteras and Ibsen (2019)[42] Volume closure, six tests Geometry, penetration, final heave, and before/after void ratio are direct; area ratio is derived. Closure uses the reported void-ratio change. Ragni et al. (2020)[15] Formula reduction, one test Prototype geometry, penetration and heave ratios, geometric decomposition, and PIV strain are direct. The retained test uses the reported self-weight boundary condition. Kim and Kim (2020)[16] Final state prediction, four profiles Geometry, final heave, suction–depth curves, and most soil properties are direct; missing state inputs follow cited analog or common rules. Digitized normalized curves are converted at prototype scale; kr/k|=1.3 29035_ 29042 68408078 29035_ 69640972 12349 28721 314 28723. Tran and Randolph (2008)[14] Final state prediction, five profiles Geometry, suction curves, and most soil properties are direct; initial state is derived where required. Final heave is calculated as hea|e=L−zend 29032_ 29032 29029 29025 69640972 29029 12349 29004 8704 29050_ 29029 29038 29028, where L 29004 is the caisson length and zend 29050_ 29029 29038 29028 is the final penetration depth digitized from the published curve. Tjelta (1995)[41] Full-scale process prediction, A3 Field suction and heave are direct; soil-state inputs combine site/design information and cited regional analogs. The solver uses ϕmodel=38.5∘ 28958_ 29037 29039 29028 29029 29036 12349 28723 28728 314 28725 8718; no finite outer boundary is imposed. Wang et al. (2019)[44] Laboratory process prediction Geometry, process curves, and the principal soil properties are direct or derived from reported values. Case inputs use the common hydraulic-ratio and prediction sensitivity band rules in S4.2. Zhao et al. (2025)[43] Model-tank process prediction, three tests Geometry, process curves, relative density, permeability, unit weight, and friction angle are direct; void-ratio bounds use a cited same-group analog. The initial void ratio is derived from the reported relative density and adopted bounds. S4.5 Benchmark Data and Summary Metrics S4.5.1 Formula Reduction and Volume Closure Formula reduction and volume closure assess whether the overall volume relation reproduces final plug heave from the end-state information reported for each case. The six Koteras and Ibsen (2019) cases use the reported void-ratio changes to close the volume balance, whereas Ragni et al. (2020) T3 uses the closed-form expression with its reported boundary condition. Supplementary Table 13. Formula reduction and volume closure summary. Case Test Measured heave (m) Predicted heave (m) Error (%) Source basis Koteras and Ibsen (2019) T1 45.00 49.53 10.07 Reported void-ratio closure. Koteras and Ibsen (2019) T2 47.00 48.57 3.33 Reported void-ratio closure. Koteras and Ibsen (2019) T3 44.00 47.60 8.19 Reported void-ratio closure. Koteras and Ibsen (2019) T4 42.00 54.40 29.54 Reported void-ratio closure. Koteras and Ibsen (2019) T5 49.00 55.40 13.06 Reported void-ratio closure. Koteras and Ibsen (2019) T9 57.00 58.45 2.55 Reported void-ratio closure. Ragni et al. (2020) T3 10.50 11.19 6.60 Closed-form evaluation with the reported Ragni boundary condition. Six of the seven calculated heaves differed from the measurements by 13.06% or less. The Ragni closed-form case had an error of 6.60%, whereas Koteras and Ibsen T4 had the largest error at 29.54%. S4.5.2 Final State Prediction Final state prediction compares measured final plug heave with the central prediction and predefined prediction sensitivity band for nine profiles. The four Kim and Kim endpoints are reported final heaves. The five Tran and Randolph endpoints are reconstructed from the final digitized penetration depth using hea|e=L−zend 29032_ 29032 29029 29025 69640972 29029 12349 29004 8704 29050_ 29029 29038 29028 under a plug-filled endpoint assumption. The prediction sensitivity band is obtained by varying eUNKNOWN/emax 29029_ 28672 68408078 29029_ over its declared sensitivity range while holding ν 28951 at 0.30. Supplementary Table 14. Final state prediction summary. Case Profile Measured final heave (cm) Central final heave (cm) Band min (cm) Band max (cm) Central error (%) Kim and Kim (2020) T1-1 S 25g 66.00 48.26 37.00 58.31 -26.87 Kim and Kim (2020) T1-2 S 50g 53.00 52.30 39.47 64.53 -1.32 Kim and Kim (2020) T3-1 S 50g 58.00 50.74 38.28 63.13 -12.52 Kim and Kim (2020) T4-3 S 25g 70.00 59.47 46.04 73.80 -15.05 Tran and Randolph (2008) TRAN-T1 50.85 60.53 48.07 70.46 19.04 Tran and Randolph (2008) TRAN-T2 45.04 41.63 34.20 47.79 -7.56 Tran and Randolph (2008) TRAN-T3 45.04 32.37 27.08 37.34 -28.13 Tran and Randolph (2008) TRAN-T4 59.93 59.67 47.34 72.31 -0.43 Tran and Randolph (2008) TRAN-T5 55.75 55.42 43.50 63.53 -0.58 Seven of the nine measured final heaves fell within the prediction sensitivity bands. The two outside-band profiles, Kim and Kim T1-1 S 25g and Tran and Randolph TRAN-T3, had central errors of −26.87% 8704 28722 28726 314 28728 28727\% and −28.13% 8704 28722 28728 314 28721 28723\%, respectively. S4.5.3 Continuous Process Prediction Continuous process prediction compares final plug heave and the measured heave–depth curve for five profiles. For all five profiles, final-point error is evaluated at the last measured depth. The field curves reported by Tjelta overlap substantially. The A3 trace was the clearest and most continuous over the installation depth and was therefore digitized for the process comparison. The complete curve is summarized by curve mean absolute percentage error (MAPE), root-mean-square error (RMSE), and normalized RMSE (NRMSE). Curve NRMSE is the curve RMSE normalized by the measured final heave. Because pointwise percentage errors can become large when measured heave is small during early penetration, curve MAPE is interpreted alongside final-point error, absolute RMSE, and NRMSE. Supplementary Table 15. Final-point comparison for continuous process prediction. Case Profile Measured final point (cm) Central final point (cm) Band min (cm) Band max (cm) Final-point error (%) Tjelta (1995) TJELTA-A3 15.43 16.70 14.45 18.98 8.23 Wang et al. (2019) full process 4.79 4.96 4.51 5.41 3.45 Zhao et al. (2025) ZHAO-T1 3.68 3.73 3.40 4.06 1.29 Zhao et al. (2025) ZHAO-T2 3.48 3.67 3.37 3.98 5.49 Zhao et al. (2025) ZHAO-T3 3.28 3.37 3.10 3.63 2.74 Supplementary Table 16. Curve-error metrics for continuous process prediction. Case Profile Curve MAPE (%) Curve RMSE (cm) Curve NRMSE (%) Tjelta (1995) TJELTA-A3 95.49 3.15 20.42 Wang et al. (2019) full process 19.24 0.52 10.80 Zhao et al. (2025) ZHAO-T1 36.73 0.16 4.22 Zhao et al. (2025) ZHAO-T2 119.25 0.14 3.92 Zhao et al. (2025) ZHAO-T3 13.88 0.14 4.25 All five measured final points fell within the prediction sensitivity bands, with final-point errors ranging from 1.29% to 8.23%. Across the complete curves, NRMSE ranged from 3.92% to 4.25% for the Zhao profiles, compared with 10.80% for Wang and 20.42% for Tjelta.