Continuous Geometry, Continuous Flow, Continuous Compression: A Numerical Component-Interaction Assessment for Fractional Clay Plasticity

1. Introduction

A constitutive model is often built the way an engineer builds a machine: design each part to perform well on its own testbench, then assemble the parts and expect the whole to behave as well as its best-tested component. Applied to clay plasticity, this assumption says that the yield surface, the flow rule, and the compression law can each be validated separately — against triaxial strength data, against dilatancy data, against oedometer data — and then combined without asking whether the combination itself changes anything. Fifty years of constitutive-modelling research has largely proceeded on this assumption, refining yield geometry, plastic flow, and compression laws as three separate problems, each owned by a different research community with its own benchmark tests [1,2]. The assumption is testable, but to the authors’ knowledge it has never actually been tested for these three components together. If a fractional-order operator acts directly on the yield surface to generate flow direction, and the hardening law that scales that operator is itself derived from the same state variable that also fixes the yield-surface size, then geometry, flow, and compression may not be separable parts at all — and fifty years of testing them one at a time would simply never have revealed it.

The three research lines were indeed developed largely in isolation. Yield geometry was refined against strength and dry-side data: Hvorslev-type and composite extensions [3,4,5,6,7,8,9] and later generalized, smooth surfaces [10] and continuous teardrop-type surfaces [11] were each validated against peak-strength envelopes, essentially independently of how flow direction or hardening was modelled. Constitutive evolution was refined against a different benchmark again — bounding-surface and subloading-surface formulations [12,13,14,15,16,17] were validated against stiffness-degradation and cyclic data, with no particular yield geometry or flow operator presumed.

Plastic flow was refined against a third, largely independent benchmark: stress-dilatancy data. Associated flow was rejected because it cannot match observed dilatancy in dense sands and overconsolidated clays [18,19,20]; classical non-associated models such as MIT-E3 and MIT-S1 [21,22,23], and more recently stress-fractional flow operators [24,25,26,27,28], were both validated primarily against dilatancy and stress-ratio data, largely independently of which yield surface or hardening law happened to be attached to them. The work of Sumelka [24] initiated the use of fractional flow rules to generate non-normality; Sun and co-workers extended this to granular soils, bounding-surface frameworks, true triaxial loading, cyclic behavior, interfaces, and overconsolidated clays; and Qu, Sun and Sumelka [27] synthesized the role of stress length scale and fractional flow in stress-fractional plasticity — an extensive literature, but one that has advanced the flow operator on its own terms.

Compression was refined against a fourth benchmark, oedometer data, and reveals the same pattern most starkly. Classical critical-state models assume the e–ln p′ relationship is linear, an assumption most realistic only for normally consolidated clays over a limited stress range [2]; work on unsaturated and highly compacted clays later showed that a single linear index cannot represent behavior consistently from low to very high stress [29]. Kaewhanam and Chaimoon (2023) found, in the context of their simplified silty sand (SSS) model, that linear functions in e–ln p are accurate at high stress but inaccurate at low stress, while curved (power-law) functions show the opposite tendency, and introduced a closed-form equation joining the two within a single expression [30]; Phonchamni et al. subsequently adapted this joining construction to the virgin compression of clays as the revised Arc Joint via Optimum Parameters (AJOP) equation, reducing settlement-prediction error by up to 85% relative to a linear index and by 10–15% relative to a purely curved index [31]. AJOP is, mechanically, the same kind of continuity fix that motivated the teardrop yield surface and the fractional flow operator — yet it has so far been used only as a stand-alone settlement equation, never embedded as the hardening modulus inside a yield-surface-and-flow-rule framework.

This is where the machine-building assumption becomes questionable rather than merely convenient. Fractional plasticity does not replace the yield surface; it replaces the operator acting on the yield surface, so its output depends on the geometry on which it acts. The hardening rate is a third ingredient that cannot be waved away as background: it scales the plastic multiplier dλ appearing in every flow-rule equation, and it governs how quickly the preconsolidation size ${p\text{'}}_c$ — and hence the point on the yield surface at which the fractional operator is evaluated — evolves along a stress path. If the compression law is itself nonlinear and stress-dependent, as AJOP asserts, then two independently state-dependent quantities, the fractional order $\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$ and the hardening modulus ${\lambda }_{AJOP}\left(p\text{'}\right)$, are driven by the same underlying state along the same stress path. There is therefore a concrete mechanical channel through which geometry, flow, and compression could interact rather than simply add up — and no existing study has tested this three-way interaction directly, because no existing study needed to: each component had already been declared successful on its own testbench.

Qu, Sun, and Sumelka’s own review [27] hints at half of this picture without completing it: it emphasizes that fractional plasticity is an efficient alternative for state-dependent nonassociativity, while separately noting the importance of improved yield-surface architecture, including reduced dry-side elastic-region surfaces for geomaterials. That observation motivates pairing geometry with flow; it says nothing about compression, which is why the present study treats compression as a third, independently motivated component rather than an afterthought bolted on for symmetry.

The objectives of this paper are therefore threefold. First, it formulates a component-based constitutive design in which yield geometry, plastic flow, and the compression/hardening law are explicitly identified as the three modified components, while constitutive evolution (mapping rules) is intentionally left unmodified so the three-way interaction can be isolated. Second, it embeds the revised AJOP compression law [31] as a nonlinear, stress-dependent hardening modulus ${\lambda }_{AJOP}\left(p\text{'}\right)$, replacing the constant compression index λ used in conventional critical-state hardening, and couples it to a continuous teardrop yield surface [11] and a state-dependent fractional stress-gradient flow rule. Third, it proposes a 2×2×2 factorial numerical assessment that separates the individual, pairwise, and three-way interaction effects of geometry, flow, and compression, so that the assumption of separability — not just the interaction hypothesis — can be tested rather than taken for granted.

2. Materials and Methods

This section formulates the constitutive framework analytically; no laboratory testing was performed for this study. Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5, Section 2.6, Section 2.7, Section 2.8 and Section 2.9 present the model formulation, and Section 2.10 presents the numerical-assessment design used to test the components’ interaction (Section 3).

2.1. Component-Based Design Concept

The proposed framework decomposes a constitutive model into interacting components. In the present study, three components are modified — yield geometry, plastic flow, and the compression/hardening law — while constitutive evolution (bounding-surface or subloading-surface mapping rules) is intentionally left in its simplest conventional form, since no mapping rule is used in either the baseline or the proposed model. This keeps the comparison in Section 3 attributable to the three components under test rather than to a fourth, unmodified mechanism.

2.2. Continuous Teardrop Yield Geometry

Let p′ be the mean effective stress, q the deviatoric stress, ${p\text{'}}_c$ the size parameter of the yield surface, and M the critical-state stress ratio. The normalized mean stress and normalized deviatoric stress are defined as

$$\bar{p}={\displaystyle \frac{p\text{'}}{{p\text{'}}_c}} , \bar{q}={\displaystyle \frac{q}{M {p\text{'}}_c}}$$

(1)

Following the continuous teardrop concept, the working yield function used in this study is

$$\bar{q}=\bar{p} · {\left(\Omega · \mathrm{l}\mathrm{n}\left({\displaystyle \frac{1}{\bar{p}}}\right)\right)}^{{\displaystyle \frac{1}{\Psi }}} , 0<\bar{p} \le 1$$

(2)

where Ψ controls the skewness/curvature of the teardrop geometry and Ω adjusts the strength scale. This is the closed-form yield function of Chatwong et al. [11] (their Equation 2, $\Omega ·ln\left(p̄\right)+{\left(\eta /M\right)}^{\Psi }=0$ with $\eta =q/p\text{'}$, rearranged into the $p̄-q̄$ form above and verified directly against the published source rather than reconstructed from memory). Ψ and Ω are taken from Chatwong et al. [11]’s own calibration for the specific clay used (their Table 1; Section 3.3). The yield function is therefore written generally as

$$f_{\mathrm{T}\mathrm{D}}\left(\bar{p}, \bar{q}; \Psi , \Omega \right)=\bar{q}-\bar{p} · {\left(\Omega · \mathrm{l}\mathrm{n}\left({\displaystyle \frac{1}{\bar{p}}}\right)\right)}^{{\displaystyle \frac{1}{\Psi }}}=0$$

(3)

2.3. Classical and Non-Associated Flow

In associated plasticity, plastic strain increments are normal to the yield surface,

$$d{\epsilon }^p=d\lambda · {\displaystyle \frac{\partial f}{\partial \sigma \text{'}}}$$

(4)

where dλ is the plastic multiplier and $\sigma \text{'}$ is the effective stress tensor. Classical non-associated plasticity instead introduces a separate plastic potential g,

$$d{\epsilon }^p=d\lambda · {\displaystyle \frac{\partial g}{\partial \sigma \text{'}}}$$

(5)

which is effective but separates stress admissibility from plastic-flow direction.

Chatwong et al. [11] themselves adopt exactly this classical non-associated route: their teardrop yield surface f is paired with a separate, simpler plastic potential g (of Original-Cam-Clay form, $ln\left(p̄\right)+\eta /M=0)$, giving a closed-form, purely linear dilatancy relation $D=M-\eta$ that is independent of the teardrop shape parameters $\Psi$ and $\Omega$. This is a legitimate and economical choice, but it means the flow direction in their model does not depend on the yield-surface shape at all — the two components are decoupled by construction. The present study deliberately does not follow that route. Section 2.4 instead applies the fractional operator directly to f, so that flow direction is generated from the same shape parameters that define the yield surface, without introducing a second, independent potential function. This is not a criticism of [11]’s choice; it is a different design decision precisely because it is the one whose consequences (Section 2.8) this paper sets out to test.

This is not a purely local design choice. Xu, Shen and Sun [32] independently develop a fractional-plasticity model for a soil–structure interface that, like the present study, combines a state-dependent fractional order with a state-dependent hardening modulus and explicitly avoids introducing an additional plastic potential; their model, however, is built on a fixed Cam-clay-type loading surface rather than a curvature-controlled yield geometry such as the teardrop surface, so it does not test whether the fractional operator’s output changes when the underlying geometry itself is varied — precisely the question Section 2.8 and Section 3 address here.

2.4. Stress-Fractional Flow Operator and Stress Length Scale

Stress-fractional plasticity modifies the flow operator instead of introducing an independent plastic potential:

$$d{\epsilon }^p = d\lambda · {}_{\mathrm{F}}{\nabla }_{\sigma \text{'}}^{\alpha } f_{\mathrm{T}\mathrm{D}}$$

(6)

where ᶠ∇^α_$\sigma \text{'}$ denotes the fractional stress-gradient operator and α is the fractional order. In principal effective stress space,

$${}_{\mathrm{F}}{\nabla }_{\sigma \text{'}}^{\alpha } f={\left[D_{\sigma \text{'}1}^{\alpha } f , D_{\sigma \text{'}2}^{\alpha } f , D_{\sigma \text{'}3}^{\alpha } f\right]}^{\mathrm{T}}$$

(7)

A stress length scale (SLS) is required to make the fractional derivative well defined. A Caputo-type stress-fractional derivative is adopted provisionally,

$$D_{{\sigma \text{'}}_i}^{\alpha } f\left({\sigma \text{'}}_i\right)=\left[{\displaystyle \frac{1}{\mathsf{\Gamma }\left(1-\alpha \right)}}\right] {\int }_{{\sigma \text{'}}_i-l_i}^{{\sigma \text{'}}_i}{\left({\sigma \text{'}}_i-\tau \right)}^{-\alpha } · \left[{\displaystyle \frac{\partial f\left(\tau \right)}{\partial \tau }}\right] d\tau , 0<\alpha \le 1$$

(8)

where $l_i>0$ is the stress length scale (SLS) in the σ′_i direction. Qu, Sun and Sumelka [27] identify three conventions for defining the SLS in the fractional plasticity literature: a past-stress (“long memory”) convention, in which l_i spans from the initial stress state to the current state; a future-reference (“critical-state”) convention, in which $l_i$ spans from the current state to a future critical state; and a combined convention using both. Sun and co-workers’ own models [25,26,28] adopt the critical-state convention, and this study follows them: $l_i$ is taken as the distance, in the relevant normalized stress coordinate, from the current state to the critical state line $\eta =M$. For both the classical and teardrop geometries used here, this critical-state point coincides exactly with the phase-transformation point ${\rho ̄}^{*}$ identified in Section 3.4 (${\rho ̄}^{*}$=0.5 for the classical ellipse, ${\rho ̄}^{*}=exp\left(-1/\Omega \right)$≈0.341 for the teardrop surface with the real $\Omega$=0.93 used in Section 3, both obtained by setting $q̄=p̄$ in the corresponding classical yield curve and the teardrop yield curve of Equation (2), respectively) — so $l_i=\left|{\rho ̄}^{*}-p̄\right|$ is fully determined by the state variables already in the model, with no additional free parameter. The numerical results in Section 3 approximate this state-dependent SLS with a fixed small step h and a finite number of terms N (Section 3.1), rather than integrating the Grünwald–Letnikov sum out to the exact distance $l_i=\left|{\rho ̄}^{*}-p̄\right|$ at every step. As implemented, the approximation departs from the critical-state convention in two respects, not one: the window has a fixed length N·h rather than the state-dependent length $l_i=\left|{\rho ̄}^{*}-p̄\right|$, and it extends in the past-stress direction (below the current $p̄$) rather than toward the critical state, so the scheme realises a fixed-length window of the past-stress type while the critical-state convention above is adopted conceptually. Replacing it with the exact state-dependent bound is a specific, well-defined numerical task for future work rather than an open conceptual question (Section 2.9). The associated flow rule is recovered as the limiting case

$$\underset{\alpha \to 1}{{}_{\mathrm{F}}{\nabla }_{\sigma \text{'}}^{\alpha }} f={\nabla }_{\sigma \text{'}} f$$

(9)

2.5. State-Dependent Fractional Order

The fractional order is linked to overconsolidation because the degree of non-associativity generally increases as the soil becomes more overconsolidated. The simplest admissible form is

$$\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)=1-\chi · \left(1-{\displaystyle \frac{1}{\mathrm{O}\mathrm{C}\mathrm{R}}}\right) , 0 \le \chi <1$$

(10)

which gives α(1) = 1 for normally consolidated clay and $\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$ < 1 for OCR > 1. A bounded exponential form may also be used:

$$\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right) = {\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}} + \left(1-{\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)·\mathrm{e}\mathrm{x}\mathrm{p}\left[-k\left(OCR-1\right)\right] , 0 < {\alpha }_{\mathrm{m}\mathrm{i}\mathrm{n}} < 1, k > 0$$

(11)

with α(1) = 1 and $\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$→αmin as OCR→∞.

2.6. Nonlinear AJOP Compression Law as the Hardening Modulus

Conventional critical-state hardening assumes a constant compression index λ, i.e., a straight virgin compression line in e–ln p′ space. The revised AJOP equation of Phonchamni et al. [31], adapted from Kaewhanam and Chaimoon [30], instead represents the e–log $\sigma \text{'}$v relationship as

$$e=\Gamma -a_c · \mathrm{l}\mathrm{o}\mathrm{g}10\left({\displaystyle \frac{{\sigma \text{'}}_v}{R}}\right)-\sqrt{\theta +{\left(a_c · \mathrm{l}\mathrm{o}\mathrm{g}10\left({\displaystyle \frac{{\sigma \text{'}}_v}{R}}\right)\right)}^2}$$

(12)

where Γ, $a_c$, R and θ are the AJOP shape parameters. The subscript on a_c is used deliberately throughout this paper to avoid confusion with the fractional order α used in Section 2.4 and Section 2.5 — the two symbols are unrelated. Because the teardrop-and-fractional-flow formulation above is written in terms of mean effective stress $p\text{'}$ rather than the one-dimensional $\sigma \text{'}v$ used in oedometer calibration, $p\text{'}\approx \sigma \text{'}v$ is adopted here as a working idealisation; a proper isotropic (or K0-consistent) recalibration of Γ,$a_c$, R and θ from triaxial or isotropic compression data is left as an open item (Section 2.9).

The instantaneous (tangent) compression modulus is obtained by differentiating Equation (12) with respect to ln $p\text{'}$. Writing x = log10(p′/R),

$${\lambda }_{AJOP}\left(p\text{'}\right)=-{\displaystyle \frac{de}{d\left(\mathrm{l}\mathrm{n} p\text{'}\right)}}=\left({\displaystyle \frac{a_c}{\mathrm{l}\mathrm{n} 10}}\right) · \left[1+{\displaystyle \frac{\left(a_c · x\right)}{\sqrt{\theta +{\left(a_c · x\right)}^2}}}\right]$$

(13)

This expression provides a useful internal consistency check. As p′→R (x→0), the bracketed term vanishes and ${\lambda }_{AJOP}$ → a_c/ln 10; as p′ becomes large (x→∞), the bracketed term tends to 1 and ${\lambda }_{AJOP}$ → $2a_c/ln10$, a constant asymptotic slope. The tangent modulus therefore rises smoothly from a lower value near the reference stress to a fixed asymptotic value at high stress — which is precisely the “curved segment joined to a linear segment” behaviour that motivated AJOP in the first place [30,31], and gives a first check that Equation (13) has been derived correctly.

The choice to use ${\lambda }_{AJOP}$ as the hardening modulus, rather than as a separate post-processing correction, follows from what a hardening modulus is defined to be. In critical-state plasticity, the modulus that appears in the hardening law is, by construction, the tangent slope of whichever compression curve is taken as the model of virgin behaviour — a constant λ is not a physical assumption in its own right, but simply the tangent slope of the linear e–ln p′ model. If AJOP is accepted as a more accurate representation of that same curve across the full stress range [30,31], consistency requires using its tangent, not the tangent of the linear approximation it was built to replace, as the hardening modulus. Equation (13) is therefore not a substitution of convenience; it is the hardening modulus implied by taking AJOP, rather than a straight line, as the compression law, obtained the same way λ itself is obtained from the linear model. The remaining question, addressed numerically in Section 3, is whether this consistent choice produces a response that differs materially from the constant-λ case once it is embedded in a flow rule.

Equation (13) also has a second, physically meaningful limit that is not captured by a constant-λ line. As $x=log10\left(p\text{'}/R\right)\to -\infty$, i.e., as $p\text{'}/R\to 0⁺$, the bracketed term in Equation (13) tends to −1, so ${\lambda }_{AJOP}$→0, while Equation (12) gives e→Γ. In other words, Γ acts as an upper bound on void ratio that AJOP approaches asymptotically at very low stress, rather than a value the compression curve passes through and continues beyond. A constant-λ line has no such bound: extrapolated to low stress it predicts an unbounded increase in e, which is not physically possible. This distinction becomes numerically important precisely in the low-stress, far-below-calibration region relevant to shallow soft-clay layers, and is examined quantitatively for three real clays in Section 3.4.

The conventional hardening law,

$$d{p\text{'}}_c=\left[{\displaystyle \frac{{p\text{'}}_c}{\left(\lambda -\kappa \right)}}\right] · d{\epsilon }_v^p (conventional, constant \lambda )$$

(14)

is therefore replaced by its AJOP-based, stress-dependent counterpart,

$$d{p\text{'}}_c=\left[{\displaystyle \frac{{p\text{'}}_c}{\left({\lambda }_{AJOP}\left({p\text{'}}_c\right)-\kappa \right)}}\right] · d{\epsilon }_v^p$$

(15)

in which ${\lambda }_{AJOP}$ is re-evaluated at the current ${p\text{'}}_c$ at every increment. Unlike the constant-λ case, the hardening modulus is now itself a function of the same stress state that also drives $\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$ in Equation (11), which is the mechanical basis for the three-way interaction hypothesis introduced in Section 2.8.

2.7. Constitutive Architecture

The complete framework consists of four sequential components:

State → Continuous yield geometry → Fractional flow → AJOP-based nonlinear hardening

(16)

First, the stress state is checked against the continuous teardrop yield surface,

$$F\left(\bar{p}, \bar{q}; \Psi , \Omega \right)=0$$

(17)

Second, the plastic flow direction is obtained from the fractional stress-gradient operator,

$$n_{\alpha }^p={}_{\mathrm{F}}{\nabla }_{\sigma \text{'}}^{\alpha } f$$

(18)

Third, the fractional order is evaluated from the overconsolidation state (in the numerical implementation of Section 3, the initial state, held fixed along each path — Section 3.1),

$$\alpha =\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$$

(19)

Finally, the hardening variable evolves according to the AJOP-based nonlinear law derived in Section 2.6,

$$d{p\text{'}}_c=\left[{\displaystyle \frac{{p\text{'}}_c}{\left({\lambda }_{AJOP}\left({p\text{'}}_c\right)-\kappa \right)}}\right] · d{\epsilon }_v^p$$

(20)

2.8. Geometry–Flow–Compression Interaction Hypothesis

The central hypothesis of this study is that geometry, flow, and compression are not separable. Since the fractional gradient acts on the yield function (Equations 6 and 7), and since the hardening law now scales the same plastic multiplier dλ that appears in that flow rule while also controlling how fast ${p\text{'}}_c$ — and hence the point on the yield surface being evaluated — evolves along a stress path, the three components are hypothesized to interact:

H_G×F×C:  the combined effect of geometry, flow and compression is not simply the sum of their individual and pairwise effects

(21)

This coupling is expected to act through two channels: a direct channel, in which ${\lambda }_{AJOP}\left({p\text{'}}_c\right)$ and $\alpha \left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$ are simultaneously evaluated at a state that also depends on which yield geometry is used to define ${p\text{'}}_c$; and an indirect, path-dependent channel, in which a nonlinear hardening rate changes the stress path itself, so that geometry and flow are subsequently evaluated at different states than they would be under conventional hardening. Section 3 is designed to test both channels numerically through a full factorial design.

2.9. Stated Approximations and Open Items

Three approximations should be stated explicitly before the numerical implementation. (i) The numerical results in Section 3 realise the stress-length-scale of Equation (8) as a fixed-length (N·h=0.2 in normalised stress), past-stress-direction, truncated Grünwald–Letnikov window rather than the exact state-dependent critical-state bound $l_i=\left|{\rho ̄}^{*}-p̄\right|$ (Section 2.4); implementing the exact bound is a specific numerical task rather than an open conceptual question. (ii) The fractional order α is evaluated at the initial overconsolidation state and held fixed along each path (Section 3.1); re-evaluating α from the evolving state at every increment is a well-defined extension left for future work. (iii) The AJOP parameters Γ, $a_c$, R, θ in Equation (12) were calibrated from one-dimensional oedometer data and have not yet been re-derived for the isotropic p′ used elsewhere in this framework.

2.10. Numerical-Assessment Design

The assessment is designed as a constitutive-component experiment rather than a conventional model-validation exercise: the aim is to determine whether the response produced by combining a stress-fractional flow operator with a nonlinear hardening law depends on the underlying yield geometry, and vice versa. Three factors are varied: yield geometry (classical reference surface or continuous teardrop surface), flow operator (classical integer-order gradient or state-dependent fractional stress-gradient), and hardening/compression law (conventional constant-λ or AJOP-based ${\lambda }_{AJOP}\left({p\text{'}}_c\right)$). This gives a 2×2×2 = 8-framework comparison.

Let R denote a response quantity of interest (e.g., predicted dilatancy, volumetric strain increment, or settlement error against laboratory/field data). The main effects are

$$\mathsf{\Delta }\mathrm{G}=R\left(\mathrm{F}1\right)-R\left(\mathrm{F}0\right) , \mathsf{\Delta }\mathrm{F}=R\left(\mathrm{F}2\right)-R\left(\mathrm{F}0\right) , \mathsf{\Delta }\mathrm{C}=R\left(\mathrm{F}3\right)-R\left(\mathrm{F}0\right)$$

(22)

the pairwise interaction indices are

$${\mathrm{I}}_{\mathrm{G}\mathrm{F}}=R\left(\mathrm{F}4\right)-R\left(\mathrm{F}1\right)-R\left(\mathrm{F}2\right)+R\left(\mathrm{F}0\right)$$

(23)

$${\mathrm{I}}_{\mathrm{G}\mathrm{C}}=R\left(\mathrm{F}5\right)-R\left(\mathrm{F}1\right)-R\left(\mathrm{F}3\right)+R\left(\mathrm{F}0\right)$$

(24)

$${\mathrm{I}}_{\mathrm{F}\mathrm{C}}=R\left(\mathrm{F}6\right)-R\left(\mathrm{F}2\right)-R\left(\mathrm{F}3\right)+R\left(\mathrm{F}0\right)$$

(25)

and the three-way interaction index, the standard factorial triple-contrast, is

$${\mathrm{I}}_{\mathrm{G}\mathrm{F}\mathrm{C}}=R\left(\mathrm{F}7\right)-R\left(\mathrm{F}6\right)-R\left(\mathrm{F}5\right)-R\left(\mathrm{F}4\right)+R\left(\mathrm{F}3\right)+R\left(\mathrm{F}2\right)+R\left(\mathrm{F}1\right)-R\left(\mathrm{F}0\right)$$

(26)

If $I_{GFC}$ is close to zero, the three components act additively even though pairwise interactions may still be present. If $I_{GFC}$ is non-zero, geometry, flow and compression jointly determine the response in a way that cannot be recovered from any subset of one- or two-component studies. A full benchmark would evaluate the eight frameworks along several stress paths — normally consolidated drained compression, overconsolidated drained compression, undrained compression, and a one-dimensional oedometer-type path that would exercise the AJOP component against the data type it was originally calibrated on [31] — but the present study implements the overconsolidated drained-shear case at constant mean effective stress as the primary test of the factorial design (Section 3.1), together with a second, deliberately approximate check under undrained conditions (Section 4.3); the remaining paths — normally consolidated drained compression and a proper oedometer-type path — are left as future work (Section 5) rather than attempted here.

3. Results

3.1. Numerical Implementation

The eight configurations of Table 2 were evaluated by incrementally integrating a drained shear path at constant mean effective stress, starting from a specified initial OCR. The shear-strain increment dεₛᵖ is treated as the controlled (independent) loading variable, and each step updates $dp̄=-p̄·D\left(p̄\right)/\left({\lambda }_{e}ff\right)·d\epsilon ₛᵖandd\epsilon ᵥᵖ=D\left(p̄\right)·d\epsilon ₛᵖ$. Two points of convention are worth stating explicitly here, since they determine how the numbers in Section 3.3 and Section 3.4 should be read. First, p′ itself is held fixed throughout a given path; it is $p̄=p\text{'}/p\text{'}c$ that evolves, and only because $p\text{'}c$ evolves under the hardening law — the path is therefore driven entirely by plastic hardening/softening at constant mean effective stress, not by an externally imposed change in $p̄$ or in $p\text{'}$ itself. Second, positive $D$ corresponds to net contraction and negative $D$ to net dilation, consistent with εᵥᵖ>0 meaning compressive volumetric strain throughout this paper (Section 2); $dp̄$and $D$ therefore carry the same sign in the update equation above, and p̄ increases when $D$ is negative (dilative) and decreases when $D$ is positive (contractive) for a fixed $p\text{'}$. The clamps $p{̄}_{min}$=0.03 and $p{̄}_{max}$=0.985 in the update equation are numerical safety bounds only — in the runs reported in Section 3.3 and Section 3.4, every framework is net-dilative from the starting state used here (Section 4.1) and so $p̄$ increases rather than decreases; each framework instead stabilizes at an intermediate asymptote well short of $p{̄}_{max}$ (specifically, its own flow rule’s critical-state point: $p̄$=0.5 for classical-integer, ${\rho ̄}^{*}$=0.341 for teardrop-integer, and the wet-side values $p̄$≈0.868 and$p̄$≈0.818 for the classical-fractional and teardrop-fractional combinations respectively — equilibria of the fixed-window Grünwald–Letnikov operator itself rather than of the underlying yield surfaces, Section 2.9 and Section 3.1.1), and none reach either clamp; $p{̄}_{min}$ and $p{̄}_{max}$ are retained as bounds for completeness rather than because it is reached by these particular paths.

where $D\left(p̄\right)$ is the flow ratio for the chosen geometry×flow combination (Section 3.1.1) and ${\lambda }_{e}ff=\lambda -\kappa or\lambda AJOP\left(p\text{'}c\right)-\kappa$ depending on the compression choice, with p′c updated each step from p̄ at fixed p′. Controlling dεₛᵖ rather than dεᵥᵖ matters because any scheme that divides by $D$ is singular wherever $D$ passes through zero: for teardrop geometry combined with integer-order (Chatwong et al. [11]’s own g-based) flow, $D\left(p̄\right)=1-{\left(\Omega ·ln\left(1/p̄\right)\right)}^{\left(1/\Psi \right)}$ passes through exactly zero at ${p̄}^{*}=exp\left(-1/\Omega \right)$ (${p̄}^{*}$≈0.341 for the real $\Omega$=0.93 used in Section 3) — a genuine phase-transformation point at which the flow direction is purely deviatoric — and dividing by $D$ there is singular. Controlling dεₛᵖ instead of dεᵥᵖ removes this singularity entirely, since $D$=0 simply means zero volumetric increment at that step rather than an undefined ratio; this is not a numerical workaround so much as the physically appropriate choice of independent variable for a strain-controlled shear test, and it is retained as the standard scheme throughout Section 3.

The response quantity is $R=\epsilon ᵥᵖ$ accumulated over a fixed shear-strain budget (εₛᵖ=6, $n_{steps}$=60,000 for full convergence, Section 3.4), so that all eight frameworks are compared after exactly the same amount of imposed shearing. This is a full incremental integration of the coupled flow rule and hardening law along a stress path, not a single-state evaluation, though it remains a simplified one-dimensional (constant-p′) idealisation rather than a general-stress-path implementation, and elastic strains before first yield are not separately tracked.

The two update equations above also imply an exact structural identity that is used throughout as an internal consistency check. Dividing them gives $d\epsilon ᵥᵖ=-{\lambda }_{eff}·d\left(lnp̄\right)$ at every increment, so for the constant-λ configurations the accumulated response has the closed form $R=-\left(\lambda -\kappa \right)·ln\left(p{̄}_{e}nd/p̄₀\right)$, where $p{̄}_{end}$ is the equilibrium point at which the configuration’s flow rule satisfies $D$=0; for the AJOP configurations, $R=-\int \left({\lambda }_{AJOP}\left(p\text{'}c\right)-\kappa \right)d\left(lnp̄\right)$ over the same interval. The flow rule therefore influences $R$ only through the location of its equilibrium point (and through whether the strain budget suffices to reach it), while the compression law influences R only through the integrand. The reported endpoints reproduce every entry of Table 4 to four to five significant figures through this identity, and — because the critical-state stress ratio M cancels from it — the results of Section 3.3, Section 3.4 and Section 3.5 are independent of M; all calculations are accordingly performed in the normalised coordinates of Equation (1) with M=1 (Section 3.1.1).

3.1.1. Flow Ratio by Geometry × Flow Combination, with Explicit Sign Derivation

For any yield function written as $f=\overline{q}-g\left(\overline{p}\right)=0$, associated flow is normal to f: the plastic strain increment is proportional to $\left(\partial f/\partial \overline{p},\partial f/\partial \overline{q}\right)=\left(-g\text{'}\left(p\right),1\right)$. The dilatancy ratio is the ratio of the volumetric to the shear component of that normal, $D=d\epsilon ᵥᵖ/d\epsilon ₛᵖ=\left(\partial f/\partial \overline{p},\partial f/\partial \overline{q}\right)=-g\text{'}\left(\overline{p}\right)$, not +g′(̄p); the minus sign is not a convention but a direct consequence of $\partial f/\partial \overline{p}=-g\text{'}\left(\overline{p}\right)$ for this form of f. This can be checked against the standard critical-state result for Modified Cam Clay written with M retained explicitly: for $f=q²-M²p\text{'}\left(p{\text{'}}_c-p\text{'}\right)=0$, the associated dilatancy ratio is $D=M²\left(2p\text{'}-p{\text{'}}_c\right)/\left(2q\right)$, which in terms of the stress ratio η=q/p′ reduces exactly to $D=\left(M²-\eta ²\right)/\left(2\eta \right)$ — negative (dilative) for η>M (the dry side, low ̄p) and positive (contractive) for η<M (the wet side), matching standard critical-state theory (e.g., Roscoe and Burland). Substituting the normalized form$̄q=g\left(\overline{p}\right)=\surd \left(\overline{p}-{\overline{p}}^2\right)$ and taking $D=-g\text{'}\left(\overline{p}\right)$ reproduces this result up to the constant factor M — in the normalised coordinates of Equation (1), $-g\text{'}\left(\overline{p}\right)=\left(M²-\eta ²\right)/\left(2M\eta \right)$ — confirming the sign: for classical (MCC) geometry with integer-order (associated) flow, $D=-g\text{'}\left(\overline{p}\right)=\left(2\overline{p}-1\right)/\left(2\surd \left(\overline{p}-\overline{p}²\right)\right)$. Consistent with this, all D values in this study are reported in these normalised units (M=1); by the closed-form identity of Section 3.1 this leaves every response in Table 4 and Table 5 unchanged. The same sign applies to the fractional-order case: the fractional operator generalizes $\partial f/\partial \overline{p}=-g\text{'}\left(\overline{p}\right)$, so D for classical geometry with fractional flow is the negative of the Grünwald–Letnikov fractional derivative of g(̄p), not its raw (positive) value. For teardrop geometry with fractional flow — the present study’s own proposal — the same reasoning applies to the real teardrop yield function$̄q=g\left(\overline{p}\right)=\overline{p}·{\left(\Omega ·ln\left(1/\overline{p}\right)\right)}^{\left(1/\Psi \right)}$: D is the negative of the Grünwald–Letnikov fractional derivative of $g\left(\overline{p}\right)$, with no separate potential required. For teardrop geometry with integer-order flow, D is instead Chatwong et al.’s [11] own $D=M-\eta$ (Section 2.3), since that is the flow rule actually proposed alongside their yield surface, and it is not derived from $-g\text{'}\left(\overline{p}\right)$ at all — it is [11]’s own non-associated potential, taken directly from their published formulation and left as published throughout this study.

Although the analytical operator in Equation (8) is written in Caputo-type form, the present numerical implementation evaluates the corresponding fractional gradient using a finite-term Grünwald–Letnikov approximation with step size h=0.005 and N=40 terms, which is commonly used for numerical fractional differentiation and converges to the same operator as the step size shrinks and the number of terms grows. The implications of replacing this approximation with the exact state-dependent integration bound now specified in Equation (8) are discussed in Section 2.9.

3.2. Compression-Law Comparison Using Real AJOP Parameters

Table 3 uses the AJOP parameters Γ, $a_c$, R and θ fitted from laboratory oedometer data for three clays by Phonchamni et al. [31] (their Table 1: Boston blue clay, London clay, Bangkok clay). For each soil, a linear compression index ${\lambda }_{LF}$ was calibrated to match AJOP exactly over R/2–2R, a typical oedometer test range straddling the reference stress — the same procedure a real laboratory calibration would follow. ${\lambda }_{AJOP}$ is reported both at p′=R (Equation 13 at x=0) and at its high-stress asymptote $2a_c/ln10$ (Section 2.6).

For all three soils, ${\lambda }_{LF}$ calibrated over R/2–2R equals ${\lambda }_{AJOP}$ evaluated at $p\text{'}=R$ exactly (Table 3), a consequence of the symmetry of AJOP’s √$\left(\theta +\left(a_{c}x\right)²\right)$ term about x=0: the secant of an even function taken symmetrically about its own center exactly recovers the function’s slope at that center. A linear index fit over a realistic test range therefore does not approximate AJOP’s high-stress asymptotic slope, which is twice as large (Table 3); it reproduces AJOP’s local tangent modulus at the reference stress R exactly, and diverges from AJOP on both sides of R. Extrapolated below the calibration range, the constant-λ line over-predicts the local compressibility relative to AJOP, and does so increasingly as p′ falls further below R, because (per the Γ-bound discussed in Section 2.6) ${\lambda }_{AJOP}$→0 as $x=log10\left(p\text{'}/R\right)\to -\infty$, while ${\lambda }_{LF}$ stays constant. This behaviour is qualitatively consistent with the settlement-error asymmetry reported by Phonchamni et al. [31] (larger linear-index error at shallow, lower-stress depths than at deep, higher-stress depths), although the present comparison is illustrative of the mechanism and is not a reproduction of their reported error magnitudes.

Of these three soils, Boston Blue Clay is used as the primary case for the factorial evaluation in Section 3.3 and Section 3.4, in preference to the other two, because it is also one of the clays for which Chatwong et al. [11] independently calibrated and validated the teardrop yield-surface parameters $\Psi$ and $\Omega$ (their Table 1). Using this soil therefore lets every parameter entering Section 3.3—the yield-surface shape, the flow rule’s implicit non-associated relation, and the compression law—be drawn from experimentally-calibrated sources rather than assumed.

As a basic sanity check on the AJOP parameters themselves, $a_c$ for Boston Blue Clay (0.430, Table 3) is itself a conventional compression index Cc in the e–log₁₀p′ sense used elsewhere in the literature, since AJOP is formulated directly in log₁₀(p′/R) terms (Section 2.6); this is a separate quantity from the ln(p′)-based λ_LF reported in Table 3, which differs from Cc by a factor of ln(10). Reported Cc values for Boston Blue Clay in the literature are commonly of the order of 0.3–0.5; the value used here falls within that range rather than outside it, which is a minimal but real check that the AJOP fit is not producing an unrealistic compression index for this specific, well-studied clay, even though it does not by itself validate the incremental framework built on top of it in Section 3.3 and Section 3.4.

Figure 1. AJOP tangent compression modulus ${\lambda }_{AJOP}\left(p\text{'}\right)$ (Equation 13) compared with a constant linear index λ calibrated over the typical oedometer test range R/2–2R, using the real Boston Blue Clay parameters of Table 3. ${\lambda }_{AJOP}$ rises smoothly from $a_c/ln10$ near the reference stress R to its high-stress asymptote 2a_c/ln10, whereas the constant-λ line stays flat across the whole stress range.

Figure 1. AJOP tangent compression modulus ${\lambda }_{AJOP}\left(p\text{'}\right)$ (Equation 13) compared with a constant linear index λ calibrated over the typical oedometer test range R/2–2R, using the real Boston Blue Clay parameters of Table 3. ${\lambda }_{AJOP}$ rises smoothly from $a_c/ln10$ near the reference stress R to its high-stress asymptote 2a_c/ln10, whereas the constant-λ line stays flat across the whole stress range.

3.3. Factorial Results (F0–F7) for Boston Blue Clay

The 2×2×2 factorial design of Table 2 was evaluated using the real Boston Blue Clay AJOP parameters from Table 3 and the real teardrop yield-surface parameters of Chatwong et al. [11] (their Table 1: $\Psi$=1.4, $\Omega$=0.93). The conventional (constant-λ) branch uses λ=${\lambda }_{LF}$=0.187, the R/2–2R secant of the same AJOP curve (Table 3) — exactly the linear index a laboratory calibration of this soil over a realistic test range would produce — and the swelling index κ=0.0365 is obtained by applying Chatwong et al.’s [11] calibrated ratio κ/λ=0.036/0.184=0.196 to this ${\lambda }_{LF}$. A representative preconsolidation state $p\text{'}c$=600 kPa and OCR=5 was used. The state-dependent fractional order (Equation 11, αmin=0.55, k=0.40) gives α(OCR=5)=0.641. Table 4 reports $R=\epsilon ᵥᵖ$ accumulated over the fixed shear-strain budget described in Section 3.1; the calculation is reproduced with live formulas in the accompanying spreadsheet (FF_calculator.xlsx).

Figure 2. The 2×2×2 factorial design of Table 2, shown as a cube. Each axis toggles one component (Geometry, Flow, Compression) from its baseline (classical / integer-order / conventional, vertex F0) to its modified state (teardrop / fractional / AJOP, vertex F7). Edges connect configurations that differ in exactly one component.

Figure 2. The 2×2×2 factorial design of Table 2, shown as a cube. Each axis toggles one component (Geometry, Flow, Compression) from its baseline (classical / integer-order / conventional, vertex F0) to its modified state (teardrop / fractional / AJOP, vertex F7). Edges connect configurations that differ in exactly one component.

Figure 3. The classical Modified Cam Clay yield curve and the real teardrop yield curve of Chatwong et al. [11] in normalized stress space, shown here for Kaolin Clay ($\Psi$=1.5, $\Omega$=1.17, Table 1 of [11]) because its parameters give the clearest teardrop shape for illustration. This figure is illustrative only; all factorial calculations in this study, including Table 4, use Boston Blue Clay’s own parameters (Ψ=1.4, Ω=0.93) throughout. The state point $\overline{p}₀$=1/OCR=0.2 (OCR=5) marked here is the same initial state used to evaluate the dilatancy ratio D in Section 3.1, independent of which soil’s curve is plotted.

Figure 3. The classical Modified Cam Clay yield curve and the real teardrop yield curve of Chatwong et al. [11] in normalized stress space, shown here for Kaolin Clay ($\Psi$=1.5, $\Omega$=1.17, Table 1 of [11]) because its parameters give the clearest teardrop shape for illustration. This figure is illustrative only; all factorial calculations in this study, including Table 4, use Boston Blue Clay’s own parameters (Ψ=1.4, Ω=0.93) throughout. The state point $\overline{p}₀$=1/OCR=0.2 (OCR=5) marked here is the same initial state used to evaluate the dilatancy ratio D in Section 3.1, independent of which soil’s curve is plotted.

Table 4. Response $R=\epsilon ᵥᵖ$ (accumulated over εₛᵖ=6) for the eight constitutive configurations (Boston Blue Clay, real parameters, OCR = 5).

Framework	Geometry	Flow	Compression	R = εᵥᵖ
F0	Classical	Integer-order	Conventional	−0.13765
F1	Teardrop	Integer-order (ref. [11])	Conventional	−0.08024
F2	Classical	Fractional	Conventional	−0.22057
F3	Classical	Integer-order	AJOP	−0.10821
F4	Teardrop	Fractional	Conventional	−0.21163
F5	Teardrop	Integer-order (ref. [11])	AJOP	−0.07502
F6	Classical	Fractional	AJOP	−0.13368
F7	Teardrop	Fractional	AJOP	−0.13190

Table 4. Response $R=\epsilon ᵥᵖ$ (accumulated over εₛᵖ=6) for the eight constitutive configurations (Boston Blue Clay, real parameters, OCR = 5).

Framework	Geometry	Flow	Compression	R = εᵥᵖ
F0	Classical	Integer-order	Conventional	−0.13765
F1	Teardrop	Integer-order (ref. [11])	Conventional	−0.08024
F2	Classical	Fractional	Conventional	−0.22057
F3	Classical	Integer-order	AJOP	−0.10821
F4	Teardrop	Fractional	Conventional	−0.21163
F5	Teardrop	Integer-order (ref. [11])	AJOP	−0.07502
F6	Classical	Fractional	AJOP	−0.13368
F7	Teardrop	Fractional	AJOP	−0.13190

All eight entries in Table 4 are negative, reflecting net dilation from the initial state$\overline{p}₀$=0.2 used throughout Section 3.3, Section 3.4 and Section 3.5: every flow rule in this study predicts dilative behaviour at this starting state (Section 4.1). As shear proceeds, $\overline{p}$ rises from 0.2 toward each framework’s own critical-state point —$\overline{p}$=0.5 for classical-integer flow, ${\rho ̄}^{*}$≈0.341 for teardrop-integer flow, and $\overline{p}$≈0.868 and 0.818 for the two fractional-flow combinations (Section 3.1) — where $D$=0 and further accumulation cease. The magnitude of each entry in Table 4 therefore reflects how far the shared starting state sits from that particular framework’s own critical point, not a transition from contractive to dilative behavior partway through the path; no entry in Table 4 represents a sign error or numerical artefact.

From Table 4, the main effects and interaction indices (Equations 22–26) are:

R(F0)=−0.13765 is itself negative (net dilation, Section 4.1), so throughout Table 5 and Figure 4 and Figure 5, “% of |R(F0)|” means each quantity is divided by the magnitude |R(F0)| with its own sign preserved in the numerator, not by the signed value R(F0) itself; dividing by the signed value would flip every percentage’s sign relative to what is reported. This convention is used because it keeps the sign of each percentage matching the sign of the underlying quantity (e.g., ΔG>0 in both raw and percentage terms), which is more directly interpretable than the alternative.

Table 5. Main effects and interaction indices, Boston Blue Clay, real parameters, OCR = 5.

Quantity	Value	% of |R(F0)|
ΔG	+0.05741	+41.7%
ΔF	−0.08292	−60.2%
ΔC	+0.02944	+21.4%
$I_{GF}$	−0.04847	−35.2%
$I_{GC}$	−0.02422	−17.6%
$I_{FC}$	+0.05745	+41.7%
$I_{GFC}$	+0.01706	+12.4%

Table 5. Main effects and interaction indices, Boston Blue Clay, real parameters, OCR = 5.

Quantity	Value	% of |R(F0)|
ΔG	+0.05741	+41.7%
ΔF	−0.08292	−60.2%
ΔC	+0.02944	+21.4%
$I_{GF}$	−0.04847	−35.2%
$I_{GC}$	−0.02422	−17.6%
$I_{FC}$	+0.05745	+41.7%
$I_{GFC}$	+0.01706	+12.4%

3.4. Sensitivity to Overconsolidation Ratio

Because a single state point cannot establish whether the interaction pattern in Table 5 is a robust feature or an artefact of one particular OCR, the same eight-framework evaluation was repeated for OCR=2–10, using the real Boston Blue Clay parameters throughout. OCR=2 ($p̄₀$=0.5) is a degenerate case for the classical baseline, since p̄=0.5 is exactly the critical-state point of the Modified Cam Clay ellipse ($D_{classical}$=0 there), giving R(F0)≈0 and ill-defined percentages; it is excluded from the range reported below. For OCR=3–10, the pattern is stable and monotonic: ΔG is positive and shrinks with OCR (+94% at OCR=3 to +24% at OCR=10), ΔF is negative and shrinks in magnitude with OCR (−115% to −37%) while remaining the single largest quantity throughout, ΔC grows substantially with OCR (+1% to +43%) rather than staying negligible, $I_{GF}$ is negative and shrinks with OCR (−75% to −21%), and $I_{FC}$ — the flow×compression interaction — is consistently large (+37% to +43% over OCR=3–8) and comparable in magnitude to the main effects. $I_{GC}$ remains smaller but not negligible (up to −19%), as does $I_{GFC}$ (up to +16% over OCR=3–8); the OCR=10 values of $I_{FC}$ and $I_{GFC}$ are undefined, for the reason given at the end of this section. Robustness to yield-surface parameter choice is established by the second real soil in Section 3.5.

3.5. Cross-Check with a Second Soil: London Clay

Because every result reported so far uses a single soil, the same OCR=2–10 sweep was repeated for London Clay, the only other soil in this study’s set for which both the teardrop yield-surface parameters (Chatwong et al. [11], Table 1: Ψ=1.1, Ω=0.95; their ratio κ/λ=0.064/0.168=0.381 applied to ${\lambda }_{L}F$=0.0586 from Table 3, giving κ=0.0223) and the AJOP compression parameters (Phonchamni et al. [31], Table 1: Γ=0.850, $a_c$=0.135, R=530, θ=0.0058) are independently calibrated and real, rather than assumed. All other settings (p′c=600 kPa, εₛᵖ budget=6) were kept identical to the Boston Blue Clay case to isolate the effect of changing soil parameters alone.

The constant-λ quantities reproduce the Boston Blue Clay pattern exactly. At OCR=5, ΔF=−60.3% is again the largest single quantity, ΔG=+39.2% is positive as for Boston Blue Clay, and $I_{GF}$=−24.5% is negative and smaller in magnitude than ΔF; across OCR=3–10, ΔF runs from −115% to −37%, ΔG from +89% to +22%, and $I_{GF}$ from −42% to −16% — the same trends, with slightly smaller magnitudes, as Boston Blue Clay. The compression-related quantities, by contrast, encounter the hardening singularity identified in Section 3.4 far earlier for this soil. Because London Clay’s AJOP transition is much sharper (θ=0.0058) and its reference stress high relative to the mean stresses used (R=530 kPa), ${\lambda }_{AJOP}$(p′c)=κ already at p′c≈190 kPa: F6 crosses this singular point from OCR=4 upward and F7 from OCR=5, so the flow×compression interaction $I_{FC}$ is defined only at OCR=3, where it is +67.2% of |R(F0)| — again the largest interaction, and larger than any Boston Blue Clay value; F3 crosses it from OCR=8 and F5 at OCR=10, so the compression main effect ΔC is defined over OCR=3–6, where it grows from +9.6% to +46.8% (compared with +1.2% to +27.8% for Boston Blue Clay over the same range), and $I_{GC}$ over the same range grows to −27%. The ranking of terms — ΔF and $I_{FC}$ largest where both are defined, $I_{GF}$ smaller and shrinking with OCR, ΔC substantial and growing with OCR — is therefore not a Boston Blue Clay-specific artefact; at the same time, the second soil demonstrates that the AJOP hardening singularity is soil-dependent and can intrude well inside the moderate-OCR range rather than only at its edge. This remains a two-soil check, not a general validation across soil types — both soils used here happen to be relatively low-plasticity, moderate-OCR marine/glacial clays for which Chatwong et al. [11] provide real parameters — and extending it to a wider range of plasticity indices remains future work (Section 5).

That the interaction pattern varies smoothly and monotonically with OCR over this range is consistent with independent microstructural evidence that OCR itself produces a graded, rather than abrupt, change in soft clay fabric: Kong et al. [33] used SEM imaging to show that the fractal dimension of soft clay pore structure decreases smoothly with increasing OCR under cyclic loading, with no discontinuity reported across the OCR range they tested. This does not validate the specific magnitudes reported here, but it is at least consistent with treating OCR as a state variable that shifts constitutive behavior gradually rather than discontinuously, which is an implicit assumption of the sensitivity analysis in this section.

This sensitivity analysis also surfaced a feature not visible in the single-state evaluation: the teardrop-plus-integer-order flow ratio $D=M-\eta$ from Chatwong et al. [11] passes through zero at ${\rho ̄}^{*}=exp\left(-1/\Omega \right)\approx$0.341 for Boston Blue Clay’s Ω=0.93, a phase-transformation point at which flow is purely deviatoric. Under the associated-flow dilatancy ratio derived in Section 3.1.1, every framework’s shear-strain-controlled path from $\overline{p}₀$=0.2 rises toward, and stabilises smoothly at, its own flow rule’s equilibrium point —$\overline{p}$=0.5 for classical-integer flow, ${\rho ̄}^{*}$ for teardrop-integer flow, and the wet-side values ̄p≈0.868 and 0.818 for the two fractional-flow combinations (Section 3.1) — and no oscillation or clamp-hugging behaviour was observed for any framework over OCR=3–8. At OCR=10, however, one framework fails structurally rather than numerically. Because ${\lambda }_{AJOP}$(p′c)→0 at low preconsolidation stress (Section 2.6), the effective hardening modulus ${\lambda }_{AJOP}$−κ necessarily vanishes at a finite p′c — for the Boston Blue Clay parameters of Table 3, at p′c≈67.4 kPa — where the hardening law of Equation (15) is singular. At OCR=10 (p′=60 kPa) this singular point corresponds to ̄p=p′/p′c≈0.891, and the classical-fractional equilibrium (̄p≈0.904 at α(10)=0.562) lies beyond it: F6’s path must cross the singularity before reaching equilibrium, and it does not converge at any step count tested — runs at $n_{steps}$=60,000–200,000 agree to three to four figures but a run at 400,000 departs from them, an apparent convergence that is spurious. The teardrop-fractional equilibrium (̄p≈0.862) lies just inside the singular point, and F7 converges normally. The compression-related interaction indices $I_{FC}$ and $I_{GFC}$, which involve F6, are therefore reported over OCR=3–8 only; every quantity not involving F6 (ΔG, ΔF, ΔC, $I_{GF}$, $I_{G}C$) is converged over the full OCR=3–10 range, and all results in Section 3.3, Section 3.4 and Section 3.5 use $n_{s}teps$=60,000. This hardening singularity is a structural property of embedding AJOP as the hardening modulus: the same low-stress flattening that gives AJOP its Γ-bound (Section 2.6) guarantees that its tangent falls below the swelling index at sufficiently low p′c. Section 3.5 shows that for a second soil it intrudes at much lower OCR, and Section 4.6 shows that a proportional-κ variant of the hardening law removes it entirely.

4. Discussion

The robust, path-integrated results of Section 3.3, Section 3.4 and Section 3.5, obtained using real yield-surface parameters from Chatwong et al. [11] and real AJOP parameters from Phonchamni et al. [31], give the following picture. The flow main effect ΔF is the largest single quantity at every OCR tested (−37% to −115% of |R(F0)| across OCR=10–3), followed by the geometry main effect ΔG (+24% to +94%, shrinking with OCR) and the flow×compression interaction $I_{FC}$ (+37% to +43% for Boston Blue Clay over OCR=3–8, roughly stable; +67% for London Clay at OCR=3, the only overconsolidation ratio at which it is defined for that soil — Section 3.4 and Section 3.5). The geometry×flow interaction $I_{GF}$ is negative throughout and shrinks in magnitude as OCR increases (−21% to −75%). The compression main effect ΔC, rather than staying small, grows substantially with OCR (+1% to +43% for Boston Blue Clay across OCR=3–10; +10% to +47% for London Clay across OCR=3–6, Section 3.5) and is comparable in magnitude to ΔG at high OCR. No single main effect or interaction dominates every other throughout the range.

This is a more interesting result than the framework’s original motivation anticipated. Section 1 motivated compression as a third component deserving the same constitutive status as geometry and flow; Section 3 confirms that expectation more strongly than the paper originally set out to claim. Compression’s main effect ΔC grows with OCR rather than staying negligible, and its interaction with flow, $I_{FC}$, is consistently one of the two or three largest quantities in the entire factorial design — comparable in magnitude to the geometry×flow interaction $I_{GF}$, and larger than it from OCR=5 upward. AJOP’s own reported strength, large settlement-error reduction at shallow depth [30,31], is a statement about compression-alone behavior under oedometer loading, a different loading path and response quantity from the shear-dominated one evaluated here; the present result is consistent with that strength rather than independent of it — compression’s influence propagates into the shear response primarily through its interaction with the flow rule ($I_{FC}$), not as a standalone effect (ΔC alone is smaller than $I_{FC}$ at every OCR tested), which is itself a specific, falsifiable claim about how AJOP’s compression law couples to the rest of this framework. The same embedding also carries a structural cost that Section 3.4 and Section 3.5 make explicit: because ${\lambda }_{AJOP}$→0 at low stress, the effective hardening modulus ${\lambda }_{AJOP}$−κ necessarily vanishes at a finite, soil-dependent preconsolidation stress (≈67 kPa for Boston Blue Clay, ≈190 kPa for London Clay), beyond which Equation (15) is singular — physically, the point at which the virgin tangent compressibility falls below the recompression index. The Γ-bound presented in Section 2.6 as AJOP’s low-stress advantage over a constant λ therefore has a precise flip side when AJOP is embedded as a hardening modulus rather than used as a settlement equation, and any sufficiently dilative path can reach it. A regularisation, or an explicit restriction of the framework’s domain of validity to p′c above the singular stress, is required before the AJOP-hardening branch of this framework can be used predictively; Section 4.6 implements one such modification — tying the swelling index to the AJOP tangent at the calibrated ratio κ/λ — and finds that it removes the singularity by construction while leaving the factorial ranking intact.

MIT-E3 [22] and MIT-S1 [23] remain a relevant point of comparison for $I_{GF}$: both are complete architectures in which yield geometry and flow were calibrated together, consistent with a non-trivial geometry×flow interaction being a real feature of coupled constitutive design rather than a modelling artefact, even though $I_{GF}$ is not the largest such interaction found here. The unified hardening (UH) model of Yao, Hou and Zhou [7] ties yield-surface size to a hardening parameter, which would predict some geometry–compression coupling; the present $I_{GC}$ is small at low OCR but grows to −19% by OCR=8–10, though UH does not use the same shear-strain-controlled response quantity and so is not a direct test of the same claim. The flow×compression interaction $I_{FC}$, for which no directly comparable coupled architecture in the literature reviewed here isolates the same pairing, is the largest interaction term found in this study and is not directly anticipated by any of the models compared against in this Discussion.

On the flow side, Sun and co-workers’ stress-fractional plasticity [25,26,28] and the Qu, Sun and Sumelka review [27] evaluate the fractional operator on a fixed, usually classical, yield surface, and so could not previously have observed $I_{GF}$, whatever its true size — the present result extends the domain over which their own stated mechanism (fractional differentiation of the yield function) has actually been tested, rather than confirming or contradicting a claim that literature made. A more direct point of contact is Chatwong et al. [11] itself: their own non-associated flow rule D=M−η (Section 2.3) was found here to pass through zero at a genuine phase-transformation point ρ̄*=exp(−1/Ω), which the fractional alternative proposed in this paper (F4, F7) does not exhibit. This is a concrete, structural difference between the two flow treatments on the same yield surface, not merely a difference in numerical convenience, and it is the kind of observation that the present factorial framework was designed to surface.

These results should still be read as a first, transparent test of the interaction hypothesis, not a validated constitutive prediction ready for design use. The teardrop yield function is now the verified closed form of Chatwong et al. [11] (Section 2.2), and the flow-integration scheme is now a genuine incremental path integration rather than a single-state algebraic proxy (Section 3.1), which closes two of the three limitations identified in earlier review of this work. The third — that all results rest on one soil and one loading path — is now only partly open. A second soil (London Clay, Section 3.5) reproduces the same ranking of terms (ΔF and $I_{FC}$ largest where defined, $I_{GF}$ smaller and shrinking with OCR, ΔC substantial and growing with OCR) with somewhat different magnitudes and with the compression-related quantities bounded above in OCR by the AJOP hardening singularity (Section 3.4 and Section 3.5), and a second, approximate loading path (undrained shear, Section 4.3) shows all eight frameworks reaching stable equilibria at the same D=0 points as the drained path, with the constant-λ quantities carrying over identically and the compression-related quantities suppressed several-fold — a path effect the drained results alone do not reveal. What remains genuinely untested is sensitivity to soil plasticity across a wider range than these two moderate-plasticity soils span, stress path types beyond drained and approximately-undrained (K0, cyclic), and the exact numerical realisation of the state-dependent stress-length-scale convention specified in Equation (8), and the effect of re-evaluating the fractional order from the evolving state along each path rather than holding it at its initial-state value (Section 2.9). What the present results do establish is that no single main effect or interaction dominates the factorial design; the flow main effect and the flow×compression interaction are consistently the largest quantities, compression’s role is substantial and grows with overconsolidation ratio rather than being negligible, and the geometry×flow interaction carries over to the approximate undrained loading path unchanged — exactly so, since within these two idealised schemes the constant-λ quantities are fixed by the same D=0 equilibrium points (Section 4.3) — while compression’s role is suppressed several-fold there, a genuine path effect driven by the elastic term dominating the undrained volumetric modulus. On this evidence the geometry×flow coupling is the more path-independent of the two interactions examined here, with the caveat that its undrained persistence is structural to the two schemes compared rather than an independent test.

4.1. Where the Classical and Teardrop Flow Rules Agree and Disagree on Dilatancy

Section 3.1.1 derives the classical associated-flow dilatancy ratio as $D=-g\text{'}\left(p̄\right)$, which reduces to the standard critical-state result $D=\left(M²-\eta ²\right)/\left(2\eta \right)$. With this sign, classical MCC predicts $D$<0 (dilative) throughout the dry side ($p̄$<0.5), matching standard critical-state theory (e.g., Roscoe and Burland) and the well-documented tendency of real heavily overconsolidated clay to dilate on shearing rather than contract [1,5,19,20]. The teardrop surface’s own non-associated flow rule $D=M-\eta$ [11] gives the same sign (dilative) throughout most of the dry side too.

What the comparison shows is a genuine, more specific disagreement: the two flow rules agree in sign almost everywhere, but not because they agree on where the critical state lies. Classical MCC’s own critical-state point is at $p̄$=0.5 (where $D_{classical}$=0 identically, independent of any soil parameter); the teardrop flow rule’s critical-state point is at ${\rho ̄}^{*}=exp\left(-1/\Omega \right)$≈0.341 for Boston Blue Clay (Section 2.4), a different, soil-dependent value. Between these two points (0.341<$p̄$<0.5), the two flow rules disagree in sign: classical MCC still predicts dilative behavior there, while the teardrop flow rule already predicts contractive behavior. Outside this narrow band, the two agree — both dilative for $p̄$<0.341, both contractive for $p̄$>0.5.

Figure 6. Dilatancy ratio D versus p̄ for the classical associated flow rule (Section 3.1.1) and the teardrop surface’s non-associated flow rule $D=M-\eta$ [11] (real Boston Blue Clay parameters, $\Psi$=1.4, $\Omega$=0.93). $D$>0 is contractive, $D$<0 is dilative. The two flow rules disagree in sign only in the shaded band between their respective critical-state points (${\rho ̄}^{*}$≈0.341 and $p̄$=0.5); outside this band they agree.

Figure 6. Dilatancy ratio D versus p̄ for the classical associated flow rule (Section 3.1.1) and the teardrop surface’s non-associated flow rule $D=M-\eta$ [11] (real Boston Blue Clay parameters, $\Psi$=1.4, $\Omega$=0.93). $D$>0 is contractive, $D$<0 is dilative. The two flow rules disagree in sign only in the shaded band between their respective critical-state points (${\rho ̄}^{*}$≈0.341 and $p̄$=0.5); outside this band they agree.

The state points explored in Section 3.3, Section 3.4 and Section 3.5 ($p̄₀$=1/OCR, OCR≥3, so $p̄₀$≤0.333) fall inside the agreement zone — both flow rules predict dilative behavior at the starting state for every OCR used in this study’s factorial design. This means the substantial numerical differences between classical- and teardrop-geometry frameworks reported in Table 4 and Table 5 (ΔG, $I_{GF}$, and so on) are not explained by a sign disagreement at the starting state; they arise from the different magnitudes and different critical-state points the two flow rules predict as the state evolves under load, which is exactly the mechanism the factorial design in this paper is built to isolate. The disagreement band identified here is a genuine, quantitatively specific difference between Chatwong et al.’s [11] flow rule and the classical associated alternative, but it is a narrower and more precise finding than a blanket sign check would suggest.

4.2. Reproducing Chatwong et al.’s Own Validation for All Four Calibrated Soils

The qualitative sign check in Section 4.1 can be extended into a direct, quantitative reproduction of the validation Chatwong et al. [11] already carried out for their own model, because both halves of that comparison are available: real undrained triaxial stress paths, compiled from the primary sources Chatwong et al. [11] present in their Figure 2 (Pestana, Whittle and Gens [34] for Boston Blue Clay, Wroth and Loudon [35] for Kaolin Clay, Gens [36] for Lower Cromer till, and Gasparre [37] for London Clay), and Chatwong et al.’s [11] own computed model-prediction curves for each tested OCR — not a static yield surface read off at one state, but the full stress-path prediction their model produces when integrated with hardening, for the same soils and the same OCR values as the test data (four to eight OCR values per soil, 20 paths, 379 points in total). Both are used here exactly as reported, with no re-fitting or re-derivation on this paper’s part.

Figure 7 plots the real test data against Chatwong et al.’s [11] own prediction curves for each soil and OCR. The fit is close throughout: the mean absolute error between predicted and measured $\overline{q}$, taken over all 379 points, is 0.049 (median 0.024), and three-quarters of all points agree to within 0.05. The one systematic exception is the OCR≈2 curve in several soils, where the mean error rises to about 0.13 — visibly greater scatter than the other OCR values, consistent with the greater scatter around the prediction line visible at OCR=2 in Chatwong et al.’s [11] own published Figure 14 for London Clay. Excluding OCR≈2, the mean absolute error over the remaining 311 points is 0.031.

This close agreement is Chatwong et al.’s [11] own validation, reproduced here rather than established anew: the yield surface, its embedded flow rule, and the hardening law that together produce this prediction curves are entirely [11]’s, calibrated and validated in their own study before this paper adopted them. Chatwong et al.’s [11] own integrated prediction curves, rather than a static yield surface evaluated at one state, are the appropriate comparison here because the critical state each curve reaches depends on OCR: for Boston Blue Clay, for example, the terminal ̄p ranges from ≈0.40 at OCR=1 to ≈0.28 at OCR=8, rather than a single fixed $\overline{p}$=${\rho ̄}^{*}$=0.341 as the flow rule’s own zero-crossing point alone would suggest. A comparison against the static yield surface would compare real data to the wrong quantity; the path-integrated predictions used here avoid that mismatch.

What this section establishes, and what it does not, should be stated precisely. It confirms that the yield surface and flow rule this paper adopts from Chatwong et al. [11] are well supported by real undrained triaxial data across four independently calibrated soils — a property of the borrowed component, demonstrated by its original authors and reproduced here for transparency, not a new result of this paper. It does not validate the fractional flow rule, the AJOP-based hardening law, the 2×2×2 factorial design, or the incremental path-integration scheme in Section 3, all of which are this paper’s own contributions and remain checked only against internal consistency (Section 3) and the qualitative sign check (Section 4.1).

4.3. A Second Loading Path: Approximate Undrained Shear

All results reported in Section 3.3, Section 3.4 and Section 3.5 use a single loading path (drained shear at constant p′, Section 3.1). As a check of whether the interaction pattern also holds under a different loading condition, an approximate undrained path was implemented for the same eight frameworks at OCR=5. Undrained conditions were enforced through the standard critical-state elastic relation (constant total volumetric strain, i.e., plastic volumetric strain compensated by elastic volumetric strain via $d\epsilon {ᵥ}^e=\kappa ·dp\text{'}/\left(\left(1+e₀\right)p\text{'}\right)$, with $e₀$=2.059 — Chatwong et al.’s [11] own reported value for Boston Blue Clay, Table 1 — held fixed as an approximation) combined with the same hardening law used elsewhere in this study. Combining these two relations collapses to a single ordinary differential equation in ̄p alone — $d\left(ln\overline{p}\right)/d\epsilon ₛᵖ=-D\left(\overline{p}\right)·\left[\left(1+e₀\right)/\kappa +1/{\lambda }_{eff}\right]$ — which avoids the compounding error of tracking p′ and $p{\text{'}}_c$ as separate multiplicatively-updated quantities.

Under the dilatancy ratios of Section 3.1.1, all eight frameworks reach a stable equilibrium under this approximate undrained path. The four integer-order-flow frameworks stabilise at the same critical-state points identified in the drained results ($\overline{p}$=0.5 for F0/F3, classical; ${\rho ̄}^{*}$≈0.341 for F1/F5, teardrop), and the four fractional-flow frameworks stabilise at exactly the same wet-side equilibria as under the drained path ($\overline{p}$≈0.868 for F2/F6, $\overline{p}$≈0.818 for F4/F7). The coincidence is exact rather than approximate: both loading schemes share the same dilatancy function $D\left(\overline{p}\right)$, and every equilibrium is a root of $D$=0 — the two paths differ only in the positive factor multiplying $D$ in the update equation, which sets how fast the equilibrium is approached but cannot move it. No framework diverges toward either the $p{̄}_{min}$} or $p{̄}_{max}$ clamp under this loading path: fractional-order and integer-order flow behave alike here, each reaching a stable equilibrium regardless of yield geometry.

Because every framework stabilises within a comparable portion of the εₛᵖ=6 budget, a direct R/percentage comparison across all eight frameworks is meaningful. For the constant-λ configurations the comparison is exact by construction: both schemes obey R=−C·ln(̄p_end/̄p₀) with a path-specific positive constant C that cancels from every percentage, so ΔG=+41.7%, ΔF=−60.2% and $I_{GF}$=−35.2% reproduce the drained percentages identically. The compression-related quantities, by contrast, are strongly suppressed under this loading path but not eliminated: ΔC=+2.6%, I_GC=−2.2%, $I_{FC}$=+9.1% and $I_{GFC}$ =+0.7%, all with the same signs as the drained results but four to eight times smaller in magnitude (drained: ΔC=+21.4%, $I_{FC}$=+41.7%), because the undrained modulus is dominated by the elastic term (1+e₀)/κ≈84, which is common to all eight frameworks and dilutes the difference between the two hardening laws. Two points of numerical bookkeeping follow from the dynamics: the constant-λ quantities are converged already at $n_{steps}$=6,000, while the compression-related quantities require $n_{steps}$≈200,000 (confirmed against 600,000), because equilibrium is reached within roughly 1% of the strain budget and the AJOP-versus-constant-λ difference accrues entirely in that short transient. The geometry×flow interaction carrying over unchanged is therefore a structural property of the two schemes compared, rather than an independent test of path-independence; the several-fold suppression of compression’s role is the genuine path effect — consistent with compression’s influence operating primarily through the volumetric hardening pathway that undrained conditions largely suppress. This is a single additional loading path, not a general demonstration, and the pore-pressure and stress-path evolution this simplified single-ODE approximation omits should be treated properly before the path-dependence suggested here is considered established; that fuller undrained assessment is left for future work (Section 5).

4.4. Decomposing the Geometry Effect: Yield-Surface Shape Versus Flow-Rule Type

F1 (teardrop, Table 2) differs from F0 (classical) in two respects simultaneously, not one: the yield-surface shape changes from the MCC ellipse to the teardrop surface, and the flow rule changes from the associated normal used for F0 to Chatwong et al.’s [11] own non-associated D=M−η, which is the flow rule actually proposed alongside their yield surface (Section 2.3) rather than the associated normal to the teardrop surface itself. The reported ΔG (Table 5) is therefore a combined geometry-and-flow-rule-type effect, not a geometry effect in isolation, and the factorial design’s orthogonality claim should be read with this in mind.

To isolate the shape effect alone, a supplementary configuration was evaluated using associated flow on the teardrop surface itself — D=−g′(̄p) for the real teardrop $g\left(\overline{p}\right)=\overline{p}·{\left(\Omega ·ln\left(1/\overline{p}\right)\right)}^{\left(1/\Psi \right)}$, the same construction used for F0’s classical $D$, applied to the teardrop surface instead of [11]’s own flow rule. Comparing this configuration against F0 isolates the shape-only effect ΔG′, holding flow-rule type (associated) fixed across both geometries. Across OCR=3–10 for Boston Blue Clay, ΔG′ is +1.3% to +5.2% of |R(F0)| — consistently about 5.5% of the reported, combined ΔG (+23.7% to +94.2% over the same range). Yield-surface shape alone, in other words, accounts for only a small fraction of the geometry main effect reported in Table 4 and Table 5; the great majority of ΔG reflects the change from associated to Chatwong et al.’s [11] non-associated flow rule that accompanies the geometry change in this study’s factorial design, not the change in surface shape itself.

This does not invalidate the factorial results in Section 3.3, Section 3.4 and Section 3.5, which correctly report the effect of moving from F0 to F1 as defined in Table 2 — a specific, stated combination of geometry and flow rule, not an abstract geometry axis. It does mean ΔG and $I_{GF}$ should be read as effects of adopting the teardrop-plus-$D=M-\eta$ combination specifically, rather than of yield-surface shape in isolation, and that a full decomposition separating shape from flow-rule type across all eight frameworks — not only F0 versus F1 — is left for future work.

4.5. Sensitivity to the Fractional Order α

The state-dependent fractional order $\alpha \left(OCR\right)={\alpha }_{min}+\left(1-{\alpha }_{min}\right)exp\left(-k\left(OCR-1\right)\right)$ uses ${\alpha }_{min}$=0.55 and k=0.40 (Section 2.5); these two constants are chosen rather than fitted to data for either soil used in this study, and no source in the fractional-plasticity literature this paper draws on [24,25,26,27,28] reports values for this specific functional form, since it is this paper’s own proposal. Because these constants set $\alpha \left(OCR=5\right)$≈0.641, the value used throughout Section 3.3 and Section 3.4, it is necessary to check whether the factorial findings depend on this specific choice or hold more generally across plausible α.

Figure 8 shows the main effects and the two largest interactions recomputed at OCR=5 with α held fixed at values from 0.3 to 0.9, spanning well beyond the range $\alpha \left(OCR\right)$ actually takes across OCR=1–10 (approximately 0.56 at OCR=10 to 0.75 at OCR=3 by the formula above, approaching 1 as OCR→1). ΔG and ΔC do not depend on α at all, since F1/F3’s own flow rules do not use the fractional operator; ΔF and $I_{FC}$, which do involve fractional flow, vary substantially with α (ΔF from −72.6% at α=0.3 to −26.9% at α=0.9; $I_{FC}$ from +52.0% to +16.7% over the same range), while $I_{GF}$ is comparatively stable (−39.5% to −33.5%). Every quantity’s sign is unchanged across the entire α=0.3–0.9 range tested. The specific ranking — ΔF largest, followed by $I_{FC}$, then ΔG — holds for α up to about 0.7–0.8. This covers the α values used at OCR≥4 but sits at its upper edge at OCR=3, where α(3)=0.75; the OCR=3 drained result itself (Section 3.4, where ΔF=−115% is the largest quantity) confirms the ranking holds there. The ranking inverts at α=0.9, where ΔG and $I_{GF}$ become the largest terms instead. The qualitative claim that no single quantity dominates the factorial design (Section 3.3 and Section 3.4) is therefore robust across the tested range; the specific claim that ΔF is the single largest quantity is robust for $\alpha \left(OCR\right)$ as actually used, but would not necessarily hold under a different, unfitted choice of α_min and k that pushed α toward 0.9 or beyond.

4.6. A Proportional-κ Variant: Removing the Hardening Singularity

The hardening singularity identified in Section 3.4 and Section 3.5 arises from the combination of a stress-dependent virgin modulus ${\lambda }_{AJOP}$(p′c) with a constant swelling index κ: because ${\lambda }_{AJOP}$→0 at low stress while κ does not, their difference must vanish at a finite preconsolidation stress. This suggests a minimal, physically motivated modification that removes the singularity at its root rather than regularising it numerically: tie the swelling index to the virgin tangent at the calibrated ratio, κ(p′c)=r·${\lambda }_{AJOP}$(p′c) with r=κ/λ taken from Chatwong et al. [11] (r=0.196 for Boston Blue Clay, r=0.381 for London Clay), applied pointwise rather than once. The effective hardening modulus becomes ${\lambda }_{eff}$=(1−r)· ${\lambda }_{AJOP}$(p′c), strictly positive everywhere in the admissible domain, and the physical requirement that recompression stiffness cannot exceed virgin stiffness is enforced at every stress rather than only at the calibration point. For the constant-λ configurations the modification changes nothing — κ=r·λ_LF is the same constant as before — so F0, F1, F2 and F4, and with them ΔG, ΔF and $I_{GF}$, are identical to Section 3.3, Section 3.4 and Section 3.5.

Under this variant every compression-related quantity is defined over the full OCR=3–10 range for both soils, including the previously non-convergent F6 cases, which now converge normally (Boston Blue Clay at OCR=10: agreement to three to four figures between $n_{steps}$=60,000 and 200,000, the path stabilising unobstructed at its fractional equilibrium $\overline{p}$≈0.904). The compression-related magnitudes shrink by roughly 15–20% relative to the constant-κ results but keep their signs and trends: at OCR=5 for Boston Blue Clay, ΔC=+17.2% (constant-κ: +21.4%), $I_{FC}$=+33.6% (+41.7%), I_GC=−14.2% (−17.6%) and $I_{GFC}$ =+10.0% (+12.4%); across OCR=3–10, ΔC grows from +1.0% to +34.8% and $I_{FC}$ stays roughly stable at +28% to +34% for Boston Blue Clay, and ΔC grows from +6.0% to +41.4% and $I_{FC}$ runs from +42% to +30% for London Clay. The qualitative conclusions of Section 3.3, Section 3.4 and Section 3.5 are unchanged: no single main effect or interaction dominates, and ΔF remains the largest single quantity at every point tested except London Clay at OCR=10, where ΔC (+41.4%) marginally overtakes it (−36.8%). The variant also has an exact structural property: because λ_eff equals (1−r) times the branch modulus uniformly across all eight configurations, the factor (1−r) cancels from every percentage in the factorial design, so the variant’s normalised results are independent of r altogether (verified numerically for r=0.10–0.30).

Two conclusions follow. First, the hardening singularity of Section 3.4 and Section 3.5 is a property of the constant-κ embedding, not of the AJOP law itself: a modification that leaves AJOP untouched and changes only the swelling-index policy removes it for both soils. Second, this variant is presented for the drained factorial only; under undrained conditions κ also enters the elastic volumetric relation (Section 4.3), and a stress-dependent κ requires a consistent treatment of the elastic law before the undrained results can be repeated, which is left for future work (Section 5).

5. Conclusions

This paper formulated a three-component constitutive architecture combining the verified closed-form teardrop yield surface of Chatwong et al. [11], a state-dependent stress-fractional flow rule, and a nonlinear hardening modulus derived analytically from the AJOP compression equation, and tested the resulting geometry–flow–compression interaction hypothesis using a 2×2×2 factorial design integrated incrementally along a shear-strain-controlled stress path. Ten findings emerge. First, the AJOP compression law was shown analytically to approach Γ as a physical upper bound on void ratio at low stress, while its tangent compression modulus ${\lambda }_{AJOP}\left(p\text{'}\right)$ (Equation 13) tends to zero, unlike a constant compression index, which is bounded in neither quantity; this distinction was numerically material using real fitted parameters for three clays (Table 3). Second, no single main effect or interaction dominates the 2×2×2 factorial design: the flow main effect ΔF is the largest single quantity at every OCR tested for Boston Blue Clay (−37% to −115% of |R(F0)| across OCR=10–3), the geometry main effect ΔG is the second largest and shrinks with OCR, and the flow×compression interaction $I_{FC}$ (+37% to +43% for Boston Blue Clay over OCR=3–8; +67% for London Clay at OCR=3, Section 3.5) is consistently comparable to the geometry×flow interaction $I_{GF}$ — exceeding it in magnitude from OCR=5 upward — which is itself negative and shrinks in magnitude as OCR increases (−21% to −75% for Boston Blue Clay). Third, and unexpectedly, compression’s role is substantial and grows with OCR rather than staying negligible: the main effect ΔC rises from about 1% to 43% of |R(F0)| across OCR=3–10 for Boston Blue Clay, and its interaction with flow ($I_{FC}$) is one of the two or three largest quantities in the entire design wherever it is defined; compression does not behave as an essentially additive component. Fourth, the sensitivity analysis surfaced a genuine phase-transformation singularity in Chatwong et al.’s own non-associated flow rule ($D=M-\eta =0$ at ${\rho ̄}^{*}=exp\left(-1/\Omega \right)$) that the fractional alternative proposed here does not exhibit — a structural, not merely numerical, difference between the two flow treatments. Fifth, a second, previously unreported structural singularity was identified in the AJOP-based hardening law itself: because ${\lambda }_{AJOP}\left({p\text{'}}_C\right)$→0 at low stress, the effective hardening modulus ${\lambda }_{AJOP}$−κ vanishes at a finite, soil-dependent preconsolidation stress (p′c≈67 kPa for Boston Blue Clay, ≈190 kPa for London Clay), where Equation (15) is singular; frameworks whose flow-rule equilibrium lies beyond this point do not converge at any step count (F6 at OCR=10 for Boston Blue Clay; F6 from OCR=4 and F7 from OCR=5 for London Clay), which bounds the overconsolidation range over which the compression-related interaction indices are defined and is the precise flip side of the Γ-bound identified in the first finding; a proportional-κ variant (Section 4.6), which ties the swelling index to the AJOP tangent at the calibrated ratio κ/λ, removes the singularity by construction, leaves the constant-λ half of the design unchanged, and reproduces the same ranking with compression-related magnitudes roughly 15–20% smaller — identifying the singularity as a property of the constant-κ embedding rather than of the AJOP law itself. Sixth, the same ranking of terms reproduces for a second, independently calibrated soil (London Clay, Section 3.5) over the range where each quantity is defined: ΔF and $I_{FC}$ remain the largest quantities, $I_{GF}$ is negative and shrinks with OCR, and ΔC grows substantially with OCR (+10% to +47% over OCR=3–6, larger in magnitude than for Boston Blue Clay over the same range). Seventh, the classical associated flow rule and Chatwong et al.’s [11] own non-associated flow rule D=M−η agree in sign almost everywhere — both dilative on the dry side, both contractive on the wet side — and disagree only in a narrow band between their two different critical-state points (ρ̄*≈0.341 and p̄=0.5); the state points explored in this study’s factorial design fall outside that band, so the numerical differences between classical- and teardrop-geometry frameworks in Table 4 and Table 5 reflect differences in magnitude and evolution, not a sign disagreement at the starting state (Section 4.1). Eighth, an approximate undrained check (Section 4.3) shows that every one of the eight frameworks reaches a stable equilibrium under this loading path, at exactly the same equilibrium points as under the drained path — an exact property, since both schemes share the same dilatancy function and every equilibrium is a root of D=0; the constant-λ percentages carry over identically for the same structural reason (ΔG=+41.7%, ΔF=−60.2%, $I_{GF}$=−35.2%), while the compression-related quantities are suppressed several-fold rather than eliminated ($I_{FC}$ from +41.7% to +9.1%, ΔC from +21.4% to +2.6%), in contrast to their substantial drained-path role (Section 3.3, Section 3.4 and Section 3.5) — a genuine path effect, driven by the elastic term (1+e₀)/κ dominating the undrained volumetric modulus. This suggests compression’s interactive role operates mainly through the drained volumetric-hardening pathway, though this rests on a single additional, simplified loading path and should not yet be generalised. Ninth, Chatwong et al.’s [11] own validation of the teardrop yield surface and its embedded flow rule is reproduced directly against 379 points from 20 real undrained triaxial stress paths across all four calibrated soils (Section 4.2), using [11]’s own computed prediction curves rather than a static yield surface: agreement is close throughout (mean absolute error 0.049 in ̄q, three-quarters of points within 0.05), with more scatter at OCR≈2 consistent with [11]’s own published results. Tenth, the geometry main effect ΔG, as reported in Table 4 and Table 5, combines a yield-surface shape change with a flow-rule type change (associated to Chatwong et al.’s [11] own non-associated $D=M-\eta$); a supplementary check holding flow-rule type fixed (Section 4.4) finds that shape alone accounts for only about 5.5% of the reported ΔG across OCR=3–10, so ΔG and $I_{GF}$ should be read as effects of the specific teardrop-plus-$D=M-\eta$ combination this study’s factorial design uses, not of yield-surface shape in isolation. The primary limitations of this version are that the main factorial evaluation still rests on one primary soil (Boston Blue Clay) and one initial-OCR sweep, with a second soil (Section 3.5) and a second, approximate loading path (Section 4.3) offered only as limited cross-checks rather than full replications of Section 3.3 and Section 3.4; the numerical scheme in Section 3 only approximates the state-dependent stress-length-scale convention specified in Equation (8) rather than implementing it exactly, and holds the fractional order fixed at its initial-state value along each path (Section 2.9 and Section 3.1); and no direct comparison of this paper’s own incremental, path-integrated predictions against digitized triaxial data has yet been made — though the yield surface and flow rule those predictions are built on have Chatwong et al.’s [11] own quantitative validation against real undrained triaxial data reproduced here for all four calibrated soils (Section 4.2). Future work should extend the factorial evaluation to additional soils and stress paths, implement the state-dependent stress-length-scale convention exactly rather than by fixed-step approximation, re-evaluate the fractional order from the evolving state along each path, extend the proportional-κ variant of Section 4.6 to a thermodynamically consistent stress-dependent elastic law and to the undrained relation, extend the approximate undrained check to a proper treatment of pore pressure and stress-path evolution, and compare predicted stress-strain and dilatancy curves against published triaxial data before the interaction hypothesis can be considered validated rather than demonstrated.