Terminal Compatibility and Admissibility of a Meridional Caputo Flow Rule in Transformed-Stress Clay Plasticity

1. Introduction

When saturated clay is sheared, its skeleton can tend to contract or dilate. A constitutive flow rule describes the irreversible volumetric deformation accompanying plastic shear. Under the ideal undrained constraint used here, total volumetric deformation is zero: the plastic tendency must be balanced by an elastic volumetric increment, which changes mean effective stress. The flow rule therefore matters even when the specimen cannot change its total volume. Critical-state clay plasticity connects this response to the loading surface, elastic compressibility and hardening [1,2]. The present question is whether replacing one flow ingredient by a fractional operator preserves those connections.

The fractional ingredient considered here is a weighted integral of slopes along a prescribed stress-space meridian. Its intended constitutive role is to supply a scalar contraction–dilation ratio. This is a phenomenological choice at continuum scale: neither a grain-scale derivation of the kernel nor a time-memory law is supplied. Stress-fractional approaches have used related ideas to construct non-associated response without an independently prescribed plastic potential [3,4,5,6,7]. Our component-interaction study combined a teardrop surface, a fixed-window fractional flow approximation and nonlinear compression [8]. The subsequent closed-form study established the teardrop terminal and evaluated a terminal-directed Caputo operator, while distinguishing the surface apex from the embedded non-associated critical state [9]. Here we examine what follows when that established scalar construction is coupled to specified incremental equations.

The first physical question is where plastic volume change ceases during continued shearing. In the published non-associated base model, the radial-mapping flow rule vanishes at the transformed stress ratio $M$ [10]. A terminal chosen from the stationary point of the teardrop surface instead places the zero according to the surface-shape parameter $\Psi$. The two states need not coincide when the second surface parameter $\Omega$ is independent. Calling either point a critical state does not establish that the same limiting stress and volume response has been retained. We derive the condition for coincidence and calculate the resulting difference in an undrained terminal stress.

A second question arises in three dimensions. Soil strength depends on the relative principal stresses, so a meridional pressure–deviator relation must be connected to compression, extension and other Lode angles. The Spatially Mobilised Plane (SMP) transformation provides such a stress representation by preserving mean stress and changing the deviatoric radius [11,12,13]. A nonlocal operator samples an interval, whereas an ordinary Jacobian describes a local change. The usual local chain rule therefore cannot be assumed for a fractional stress derivative [4,14]. We identify a sufficient condition under which the present scalar integral is unchanged by the transformation, and keep that statement separate from the strain mapping required for a complete physical model.

The third question is how the selected flow changes the incremental response. We first eliminate volumetric deformation and hardening to derive the normally consolidated undrained stress locus. This shows why a close stress-path fit may carry no independent information about the flow direction. We then test two direct substitutions in the inherited overconsolidated loading structure: retaining its loading modulus, or coupling that modulus directly to the fractional dilatancy. The purpose is to expose the consequences of these precise choices before calibration. The contributions beyond [9] are the explicit terminal-compatibility condition, the exact normally consolidated locus identity, the invariant-coordinate composition criterion and the incremental substitution diagnostics. Together they specify requirements for subsequent model development; they do not constitute a validated three-dimensional fractional clay model.

1.1. Relation to Existing Models and Organization of the Study

Fractional soil models use different mathematical objects. A component-wise fractional gradient requires a multiaxial definition [5,6,7], while state-dependent scalar stress–dilatancy relations specify one part of the flow [5,6,15]. Recent formulations include shift stress, cyclic interfaces, overconsolidation, reduced elastic domains and dissipation analysis [16,17,18,19,20]; anisotropic clay models additionally introduce rotational hardening [7,21]. Bounding-surface response, non-coaxiality and calibration are separate parts of the constitutive architecture [22,23,24]. Our one-dimensional meridional integral provides none of those components by itself. Table 1 identifies the question answered by each part of this study and the evidence needed to connect it to physical prediction.

Section 2 introduces the physical assumptions, stress variables and numerical update. Section 3 asks whether the prescribed zero represents the same terminal state as the base model. Section 4 treats the stress transformation. Section 5 examines the common normally consolidated locus and the two overconsolidated substitutions. Section 6 interprets what the calculations establish about model construction and what remains to be tested experimentally.

Table 1 is the logical map of the study. Coordinate composition asks whether the same scalar operation is obtained before and after changing the stress representation. Terminal compatibility asks whether the zero selected by that operation is the state at which the inherited material model permits continuing shear without further plastic volume change. Incremental consistency asks whether the stress increment, plastic volume change and evolution of the surface size satisfy one common set of equations. The work-sign test then asks a separate question about the declared stress–plastic-coordinate pairing, whereas experimental validation requires measured observables that are actually sensitive to the chosen flow. Passing one test therefore supplies no automatic evidence for the next. This separation is used throughout the paper to prevent a mathematical property of one scalar quantity from being read as validation of a complete constitutive model.

2. Physical Basis and Mathematical Construction

2.1. Physical Assumptions and the Role of Each Equation

We work at the scale of a homogeneous material point, using an effective-stress, small-increment elastoplastic description. Elastic volumetric deformation is recoverable; plastic volumetric deformation changes the hardening state. For the ideal undrained calculations, their increments sum to zero. This provides a simple physical reading of the pressure update: with compression positive and bulk modulus $K>0$, a contractive plastic increment $d{\epsilon }_v^p>0$ requires elastic expansion and hence $dp\prime =-K\hspace{0.17em}d{\epsilon }_v^p<0$. A dilative plastic increment has the opposite effect. These implications follow from the stipulated volume balance. No pore-pressure field or drainage problem is solved.

The surface describes the stress states associated with the current hardening scale. The flow rule selects the relative plastic volumetric and shear increments. The plastic multiplier sets their magnitude, and hardening describes how the surface changes after plastic volume change. These ingredients must agree: using one volume increment to change stress and a different one implicitly to change the surface breaks the intended consistency relation. The overconsolidated tests below examine precisely that issue.

In a physical meridional model, dilatancy would be a ratio of plastic volumetric to plastic shear strain increments. Here we first define the analogous scalar-coordinate ratio $D=d{\epsilon }_v^p/d{\chi }^p$ during active plastic loading. Positive $D$ means contraction and negative $D$ means dilation within this declared coordinate system. The variable $\chi$ is introduced through a transformed elastic update and has not been proved equal to measured shear strain. A zero of $D$ is therefore called an operator zero or candidate terminal until the remaining stress, hardening and kinematic conditions are specified.

2.2. Teardrop Surface and the Two Candidate Terminals

Compression is positive. Mean effective stress $p\prime$ measures the isotropic part of stress carried by the soil skeleton. A tilde distinguishes the transformed deviatoric magnitude $\tilde{q}$ from the real magnitude $q$; all extension deviators are reported as non-negative magnitudes. Define $\tilde{\eta }=\tilde{q}/p\prime$. The current surface size $P_0$ has stress units and changes with plastic volumetric deformation. The positive dimensionless parameters $M$, $\Psi$ and $\Omega$ specify the transformed critical stress ratio and the teardrop geometry. The inherited loading surface [10] is

$$F\left(p\prime ,\tilde{q},P_0\right)={\left({\displaystyle \frac{\tilde{\eta }}{M}}\right)}^{\Psi }-\Omega \ln \left({\displaystyle \frac{P_0}{p\prime }}\right)=0.$$

(1)

At isotropic loading, Equation (1) gives the intercept $p\prime =P_0$. Thus $P_0$ records the pressure scale of the current surface, rather than a fixed normalization. With $\bar{p}=p\prime /P_0$, $\bar{q}=\tilde{q}/\left(M{P}_0\right)$ and $t=\mathrm{l}\mathrm{n}\left(1/\bar{p}\right)$, its positive meridian becomes

$$\overline{q}=g\left(t\right)=e^{-t}{\left(\Omega t\right)}^{1/\Psi }, {\displaystyle \frac{\tilde{\eta }}{M}}={\left(\Omega t\right)}^{1/\Psi }.$$

(2)

The coordinate $t$ increases as pressure decreases relative to the current surface size. It is dimensionless and is not time. In Equation (2), $\Psi$ controls the exponent and hence the curvature; $\Omega$ also scales the stress ratio at a given $t$. These are phenomenological surface parameters, not separately measured microstructural quantities. Along the surface, the base-model condition $\tilde{\eta }=M$ gives $t_b=1/\Omega$. Differentiating the surface ordinate $g\left(t\right)$ instead gives a stationary point at $t_{*}=1/\Psi$. Subscripts $b$ and $\ast$ distinguish the material terminal selected by the base flow from the surface stationary point.

The same coordinate needs a different reading at an overconsolidated initial state. Under the normalization used in the scalar diagnostics, the isotropic starting value is $t_0=\mathrm{l}\mathrm{n}\left(\mathrm{OCR}\right)$; this is a geometric statement about the position of the current stress relative to the outer-surface size. It is not a record of how long the specimen has been loaded. Likewise, the prescribed value $t_{*}=1/\Psi$ is the destination of the terminal-directed integration interval, not a reversal point in a chronological loading path. A cyclic-memory interpretation would require additional history variables and update rules that are absent here. To keep the two endpoint descriptions visible, the labels ‘op’ and ‘*’ both denote the operator-defined terminal, whereas ‘base’ and ‘b’ denote the independently selected base-model terminal.

2.3. What the Terminal-Directed Caputo Operator Represents

An ordinary associated flow rule uses the local normal to a surface. The present construction replaces one local slope by an integral of meridional slopes between the current stress coordinate and a prescribed terminal. For $0<\alpha <1$ and $t>0$, we define the candidate $D={\mathcal{D}}_{t_{*}}^{\alpha }g$ by

$${\mathcal{D}}_{t_{*}}^{\alpha }g\left(t\right)={\displaystyle \frac{1}{\Gamma \left(1-\alpha \right)}}\left\{\begin{array}{ll}{\int }_t^{t_{*}}g\prime \left(s\right){\left(s-t\right)}^{-\alpha }\hspace{0.17em}ds,& t<t_{*},\\ 0,& t=t_{*},\\ {\int }_{t_{*}}^{t}g\prime \left(s\right){\left(t-s\right)}^{-\alpha }\hspace{0.17em}ds,& t>t_{*},\end{array}\right. t_{*}={\displaystyle \frac{1}{\Psi }}.$$

(3)

The kernel weights slopes closer to the evaluation point more strongly, with $\alpha$ controlling that weighting. The chosen orientation gives a positive value on the $t<t_{*}$ side and a negative value on the $t>t_{*}$ side. The interval collapses at $t_{*}$, where the value is prescribed as zero. These are construction properties. Interpreting the output as a volumetric-to-shear plastic-coordinate ratio is the additional constitutive assumption tested below. During each evaluation, $P_0$, $\Psi$, $\Omega$ and $M$ are held fixed. The sampled points belong to the current meridian; they need not be states previously visited by the specimen. In particular, Equation (3) contains no elapsed time, loading rate or stored reversal history.

Standard fractional definitions are given in [25,26], and the earlier closed-form construction is reported in [9]. Under its notation, the present sign convention is ${\mathcal{D}}_{t_{*}}^{\alpha }g=-D_{\mathrm{c}\mathrm{a}\mathrm{p}}$; the earlier pressure-coordinate evolution contains the additional positive factor $\mathrm{e}\mathrm{x}\mathrm{p}\left(t\right)$, and recovering physical dilatancy introduces $M$. Those factors preserve the zero and sign of the present integral but change its magnitude. We deliberately test Equation (3) directly as $D$, so our trajectories must not be read as reproductions of that earlier evolution law.

As $\alpha \to 1^{-}$, Equation (3) tends to $g\prime \left(t\right)$ at an interior point. The associated pressure-coordinate slope would instead involve $M\mathrm{e}\mathrm{x}\mathrm{p}\left(t\right)g\prime \left(t\right)$. Thus the integer-order limit of this integral is well defined, but does not by itself recover the complete associated constitutive model. The reference value $\alpha =0.641$ is inherited, not calibrated against the present records. At $t=0$, the integral is finite only if $\alpha <1/\Psi$; all sign proofs use $t>0$, and the numerical cutoff is specified in Section 2.5.

This distinction also clarifies the limited meaning of nonlocality in the present paper. The ordinary surface normal uses the slope at the current stress coordinate, whereas Equation (3) combines slopes from the current coordinate to the prescribed terminal on the same current meridian. The parameter $\alpha$ changes the relative geometric weight assigned within that interval. Because the interval is in stress space, $\alpha$ is not interpreted here as viscosity, a rate parameter or a measure of temporal memory. Its physical identification would require observations and a complete stress–strain relation that are not part of this audit.

2.4. Why Transformed Stress Is Introduced

For positive principal effective stresses ${\sigma }_i^{\prime }$, let $I_1=\sum _i{\sigma }_i^{\prime }$, $I_2=\sum _{i\text{<}j}{\sigma }_i^{\prime }{\sigma }_j^{\prime }$ and $I_3={\sigma }_1^{\prime }{\sigma }_2^{\prime }{\sigma }_3^{\prime }$, and define $T_1$ = $I_1{I}_2$ − $I_3$ and $T_9$ = $I_1{I}_2$ − 9$I_3$. Write $\mathit{s}={\mathit{\sigma }}^{\prime }-p^{\prime }\mathit{I}$ and $q=\sqrt{3\mathit{s}:\mathit{s}/2}$, where $\mathit{I}$ is the identity tensor. The operational SMP-equivalent deviator used in the supplied implementation, away from the hydrostatic state, is [11,12]

$$\tilde{q}={\displaystyle \frac{2{\mathbf{I}}_1}{3\sqrt{T_1/T_9}-1}}, \tilde{\mathsf{\sigma }}\prime =p\prime \mathbf{I}+{\displaystyle \frac{\tilde{q}}{q}}\mathbf{s}, \tilde{p}\prime =p\prime .$$

(4)

Equation (4) changes the deviatoric radius while preserving mean effective stress. Under triaxial compression the implemented radius equals the real deviator; under extension the inverse is $q=3p\prime \tilde{q}/\left(3p\prime +\tilde{q}\right)$. This allows a common transformed stress ratio to correspond to different real compression and extension strengths. At the hydrostatic state, both deviators vanish and the mapped stress is taken as $p\prime \mathbf{I}$. The scalar composition proof uses only the invariance of $p\prime$ and, for fixed $P_0$, of $t$. It does not require or establish differentiability of the complete tensor mapping at every hydrostatic approach.

2.5. Data Provenance and Numerical Protocol

The numerical sets in Table 2 were carried forward from the base-model calibration files associated with [10]; no parameter is fitted in the present study. The Fujinomori stress points used here were digitized from Matsuoka et al. [12], which attributes the original test to Nakai and Matsuoka [27]. The kaolin parameter lineage is the Banerjee--Stipho test family [28,29], using plots reproduced by Xu et al. [30]. Eastern Osaka is represented only by an inherited parameter set associated with the natural-clay study [31], not by a validated trajectory here. S1 identifies the files actually used and preserves additional transcriptions without treating them as evidence for the fractional candidate. In particular, the kaolin pair $\left(\Psi ,\Omega \right)=\left(1.20,0.91\right)$ is distinct from the Wroth--Loudon pair $\left(1.50,1.17\right)$ used in [8].

We evaluate Equation (3) by SciPy algebraically weighted quadrature [32]. A table with 1800 distance-coordinate points on each side of the terminal includes the exact zero and 500 extra geometric points near the isotropic cutoff. Comparison with direct quadrature at 30 seeded interior pressures and 70 pressures between 0.985 and 0.9995 gives a maximum absolute discrepancy of $1.75\times {10}^{-11}$. The evaluation is capped at $\bar{p}=0.9995$, and the table spans $0.02\le \bar{p}\le 0.9995$. This cutoff defines the numerical initial neighborhood; it is not an independently measured material parameter.

The imposed variable $\chi$ enters the scalar elastic update $d\tilde{q}=3G\left(d\chi -d{\chi }^p\right)$. We use $K=p\prime /{\kappa }_c$, $G=3K\left(1-2\nu \right)/\left[2\left(1+\nu \right)\right]$ and $d\mathrm{l}\mathrm{n}{P}_0=d{\epsilon }_v^p/\left({\lambda }_c-{\kappa }_c\right)$. Fujinomori runs use 240,000 increments to $\chi =0.15$; kaolin F1 diagnostics use 60,000 increments to $\chi =0.30$. Initial conditions are $P_0=1$, $p{\prime }_0=1/\mathrm{O}\mathrm{C}\mathrm{R}$, ${\tilde{q}}_0={10}^{-9}$ and $\nu =0.30$. Appendix B gives the multiplier, image-point evaluation, explicit update and stopping rules so that the link between flow and output is fully specified.

The incremental cause-and-effect sequence is therefore as follows. An imposed $d\chi$ supplies a transformed shear-loading increment. Together with the elastic moduli, surface-normal components, loading modulus and selected flow ratio $D$, it determines the plastic multiplier $d\Lambda$.

The resulting multiplier sets the plastic shear and volumetric-coordinate increments. The volumetric increment changes $P_0$ through hardening, while the zero-total-volume condition produces the accompanying change in $p\prime$; the transformed deviator is then updated and the next values of $t$ and $D$ are evaluated from the new state. A change in $D$ consequently propagates through pressure, hardening and the next surface location. This is why replacing the flow factor without rederiving the multiplier and hardening terms is not a local cosmetic substitution.

To compare with the experiment, calculated stresses are converted back to real triaxial extension coordinates. Measured strain is not overlaid on $\chi$: the differential of the nonlinear stress map, discussed in Section 4.3, has not been coupled here to a physical strain transformation. Two additional diagnostics are the surface residual and a declared plastic-work increment

$$r_F=\Omega \ln \left({\displaystyle \frac{p\prime }{P_0}}\right)+{\left({\displaystyle \frac{\tilde{\eta }}{M}}\right)}^{\Psi }, d{W}^p=p\prime \hspace{0.17em}d{\epsilon }_v^p+\tilde{q}\hspace{0.17em}d{\chi }^p.$$

(5)

The tests impose $d{\epsilon }_v=0$. In Equation (5), $r_F=0$ is on the surface and $r_F>0$ is outside it. The second expression is the work assigned to the declared transformed stress–plastic-coordinate pairing, evaluated at the pre-increment stress. We monitor its sign, the loading-modulus sign and the algorithmic denominator separately. A negative value fails the chosen work condition; calling that quantity physical dissipation would additionally require a strain mapping and a stored-energy model. S1 documents the reproduction environment: Python 3.12.14, NumPy 2.3.5, SciPy 1.17.0 and Matplotlib 3.10.8.

3. Does the Operator Preserve the Intended Critical State?

3.1. Exact Condition

At a critical state, continued shearing should be possible without continuing change in stress or volume within the adopted idealization. A zero of volumetric flow is therefore necessary, but by itself is not a complete critical-state condition. We first locate the zero imposed by the scalar construction. Differentiating Equation (2) gives

$$g\prime \left(t\right)={\Omega }^{1/\Psi }{e}^{-t}{t}^{1/\Psi -1}\left({\displaystyle \frac{1}{\Psi }}-t\right).$$

(6)

The derivative changes sign at $t_{*}=1/\Psi$, which is the zero anchoring Equation (3). On the positive meridian, its stress coordinates are

$${\displaystyle \frac{{\tilde{\eta }}_{\mathrm{op}}}{M}}={\left({\displaystyle \frac{\Omega }{\Psi }}\right)}^{1/\Psi }, {\overline{p}}_{\mathrm{op}}=e^{-1/\Psi }.$$

(7)

In the base radial-mapping model, the normally consolidated rule is $D_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}=M-\tilde{\eta }$. Continued plastic shear therefore has zero volumetric increment at $\tilde{\eta }=M$, which gives ${\bar{p}}_b=\mathrm{e}\mathrm{x}\mathrm{p}\left(-1/\Omega \right)$. Comparing that independently selected state with Equation (7) yields

$$\left({\tilde{\eta }}_{\mathrm{op}},{\overline{p}}_{\mathrm{op}}\right)=\left({\tilde{\eta }}_{\mathrm{base}},{\overline{p}}_{\mathrm{base}}\right) \iff \Omega =\Psi .$$

(8)

The two definitions select the same point if and only if $\Omega =\Psi$, independently of fractional order. The physical issue is that the zero specifies where the chosen contraction–dilation mechanism ceases. Moving it changes the candidate terminal stress, even though the yield-surface equation itself has not been changed. Because the base model is non-associated, its zero need not lie at the surface apex. Equality of the terminals also does not imply equality of the complete flow direction or of its magnitude away from that point.

3.2. Consequence for Normally Consolidated Undrained Terminals

The consequence for an undrained specimen follows from the declared volume balance. Elastic compression contributes $d{\epsilon }_v^e={\kappa }_c\hspace{0.17em}d\mathrm{l}\mathrm{n}p\prime$, while isotropic hardening gives $d{\epsilon }_v^p=\left({\lambda }_c-{\kappa }_c\right)d\mathrm{l}\mathrm{n}{P}_0$. Setting their sum to zero links a change in surface size to a change in effective pressure. For $0<{\kappa }_c<{\lambda }_c$, integration from the isotropic normally consolidated pressure $p_i$ to either candidate zero gives

$${\displaystyle \frac{p_{\mathrm{base}}}{p_i}}=\exp \left[-{\displaystyle \frac{{\lambda }_c-{\kappa }_c}{{\lambda }_c\Omega }}\right],{\displaystyle \frac{p_{\mathrm{op}}}{p_i}}=\exp \left[-{\displaystyle \frac{{\lambda }_c-{\kappa }_c}{{\lambda }_c\Psi }}\right].$$

(9)

Equation (9) evaluates the two candidate endpoints on the normally consolidated locus; reaching them requires an admissible incremental branch. The calculation depends on the ratio of the two compression coefficients, so a common specific-volume normalization cancels. The complete derivation of the locus is given in Section 5.1. Table 3 distinguishes the surface coordinate $\bar{p}=p\prime /P_0$ from the terminal mean stress normalized by the initial pressure $p_i$.

For the compatible Fujinomori pair, both constructions select the same terminal. For Eastern Osaka, their normally consolidated mean stresses differ by about 1.6%. The kaolin separation is larger: the operator zero is at $\tilde{\eta }=0.794M$, and its candidate terminal pressure is 18.6% above that selected by the base flow. Thus a change presented as only a replacement of the flow can also change the limiting undrained effective-stress state. The 18.6% is a difference between two model-defined endpoints under identical inherited coefficients; it is not a measured prediction error. Figure 1 illustrates the geometry behind this difference.

This pressure difference has a conditional undrained interpretation: if two material-point calculations share the same prescribed total mean stress, the relation $p\prime =p-u$ would assign the higher effective pressure to a lower pore pressure $u$ by the same amount. No total-stress path or pore-pressure prediction is supplied here, so Table 3 reports effective-stress endpoints only. The comparison demonstrates why the location of a volumetric-flow zero matters physically without claiming a new pore-pressure validation.

4. Can the Scalar Operator Be Carried Through a Stress Transformation?

4.1. Why the Ordinary Local-Chain Substitute Fails

A nonlinear coordinate map changes distances throughout the interval sampled by a fractional operator. Its Jacobian at the current state therefore need not reproduce the transformed weighting of all the sampled slopes. To show this directly, let $T$ be smooth and strictly increasing, let $x>a$, and let $f$ be continuously differentiable on the transformed interval. For a left-sided Caputo derivative with terminal $a$, differentiating the composite inside the integral gives

$$D_{C,a}^{\alpha }\left(f\circ T\right)\left(x\right)={\displaystyle \frac{1}{\Gamma \left(1-\alpha \right)}}{\int }_a^{x}f\prime \left[T\left(s\right)\right]T\prime \left(s\right){\left(x-s\right)}^{-\alpha }\hspace{0.17em}ds.$$

(10)

The common local-chain substitute multiplies a derivative evaluated in the transformed coordinate by the current Jacobian,

$$T\prime \left(x\right)\hspace{0.17em}{D}_{C,T\left(a\right)}^{\alpha }f\left[T\left(x\right)\right].$$

(11)

The exact expression contains the local derivative $T\prime \left(s\right)$ at every integration point and the kernel measured in the original coordinate. The substitute evaluates a different integral in the transformed coordinate and multiplies it by $T\prime \left(x\right)$ only at the endpoint. For the illustrative map $T\left(x\right)=x+0.25x^2$, $f=g\left(\cdot ;\Psi =1.4,\Omega =0.93\right)$, $a=0.6$ and $x=1.2$, independent weighted quadrature gives a relative discrepancy of 19.704% at $\alpha =0.641$ (Figure 2). The map coefficient and parameter values are illustrative choices, not values fitted to any clay; the marked fractional order is the inherited reference value stated in Section 2.3. This is a counterexample to the proposed local-chain substitute, not a 19.704% error of the SMP mapping or of a clay prediction. The discrepancy tends to zero as $\alpha$ approaches one, when the ordinary derivative is recovered.

4.2. A Sufficient Invariant-Coordinate Criterion

The failure of Equation (11) does not prevent every transformed-stress use. A narrower construction avoids a chain rule altogether. Let a scalar state function be $h\left(\mathit{\sigma })=\phi \right[\xi (\mathit{\sigma }\left)\right]$, where $\phi$ is absolutely continuous on the integration interval. Let ${\mathcal{D}}_{\xi }$ be a one-dimensional nonlocal operator acting only in $\xi$, with all other state variables held fixed, and let the stress transformation $S$ satisfy $\xi \left(S\mathit{\sigma })=\xi \right(\mathit{\sigma })$. If the integration terminal is expressed in the same invariant coordinate, then

$$\left[{\mathcal{D}}_{\xi }\right(h\text{∘}S\left)\right](\mathit{\sigma })=[({\mathcal{D}}_{\xi }h\left)\text{∘}S\left]\right(\mathit{\sigma })=\Phi \right[\xi (\mathit{\sigma }\left)\right].$$

(12)

Both sides integrate the same function over the same interval with the same kernel, because the transformation leaves the integration coordinate and terminal unchanged. This proves a sufficient criterion for the scalar construction. It is not a general chain rule for component-wise fractional gradients. In physical terms, the stress representation may change the deviatoric radius while leaving the pressure coordinate sampled by this particular operator untouched.

4.3. SMP Corollary and Its Boundary

For Equation (3), the integration coordinate is $\xi =t=\mathrm{l}\mathrm{n}\left(P_0/p\prime \right)$. Equation (4) preserves $p\prime$, and $P_0$ is held fixed as an internal scalar during the integral evaluation. Consequently, both $t$ and $t_{*}$ remain unchanged: applying the stress transformation before or after the scalar operation gives the same value. This establishes coordinate compatibility for the chosen meridional integral.

A corresponding physical strain transformation is a separate requirement. For example, differentiating the extension inverse gives $dq=\left[3{\tilde{q}}^2/{\left(3p\prime +\tilde{q}\right)}^2\right]dp\prime +\left[9p{\prime }^2/{\left(3p\prime +\tilde{q}\right)}^2\right]d\tilde{q}$. Both stress increments contribute. It is therefore not justified to reinterpret the coordinate in $d\tilde{q}=3G\left(d\chi -d{\chi }^p\right)$ as measured shear strain solely because the scalar integral is invariant. The real-stress elasticity and transformation derivatives in [12] illustrate why this extra step matters. In this study, $\chi$ remains a declared numerical coordinate.

The distinction can also be expressed through power. A stress transformation becomes constitutively complete only after a dual strain transformation is chosen so that the stress–strain power is preserved. Mean-stress invariance and equality of the scalar Caputo integral establish neither that dual relation nor work conjugacy. The exact scalar composition proved here is therefore a statement about representation of the stress-coordinate operation; it is not, by itself, a proof that $d\chi$ is a physical shear-strain increment or that the transformed stress product equals material dissipation.

5. What the Incremental Calculations Establish

5.1. Normally Consolidated Fujinomori Stress Locus and Loading Coordinate

The Fujinomori pair satisfies $\Omega =\Psi$, so terminal mismatch does not obscure the comparison. We keep the same surface, elastic coefficients, hardening, stress transformation and initial state. One run uses the inherited radial-mapping dilatancy; the other uses Equation (3) with the loading modulus coupled consistently at normal consolidation. The calculation isolates the consequences of that scalar substitution. Neither run is presented as a complete reproduction of all ingredients in a published three-dimensional model.

The effective-stress locus can be determined without solving for either plastic multiplier. Under zero total volume change, elastic volume change and plastic hardening cancel. Since $t=\mathrm{l}\mathrm{n}{P}_0-\mathrm{l}\mathrm{n}p\prime$, their differential relation becomes

$${\kappa }_c\mathrm{dln}{p}^{\prime }+({\lambda }_c-{\kappa }_c)\mathrm{dln}{P}_0=0,\mathrm{dln}{P}_0=\mathrm{d}t+\mathrm{dln}{p}^{\prime }\text{⇒}\mathrm{dln}{p}^{\prime }=-{\displaystyle \frac{{\lambda }_c-{\kappa }_c}{{\lambda }_c}}\mathrm{d}t.$$

(13)

This equation says that pressure loss and growth of the surface size are linked by the prescribed compression coefficients. Substitution into the surface relation and integration from $p\prime =P_0=p_i$ give

$${\displaystyle \frac{p^{\prime }}{p_i}}=e^{-at},{\displaystyle \frac{\tilde{q}}{p_i}}=M{e}^{-at}(\Omega t{)}^{1/\Psi },a={\displaystyle \frac{{\lambda }_c-{\kappa }_c}{{\lambda }_c}}.$$

(14)

No flow ratio appears in Equation (14). Under these assumptions, changing the flow ratio changes the progression along an admissible reachable branch, while the equations fixing that branch in stress space remain the same. Accordingly, Table 4 gives nearly identical effective-stress-path RMS discrepancies, 4.31200% and 4.31311% of $p_i$. The difference of about 0.0011 percentage point reflects numerical integration and curve sampling along their common locus. The shared agreement in Figure 3a therefore tests the adopted stress locus; it cannot identify the fractional order or validate one flow candidate over the other.

The progression coordinate does change. The first crossing of 99% of the analytical terminal occurs at $100\chi \approx 1.45$ for the base candidate and 2.99 for Equation (3), as shown in Figure 3b. Appendix A explains this through two contributions: the elastic shear increment and the plastic shear coordinate required to balance volumetric hardening. The target is the analytical terminal, not the deviator reached at the imposed 15% endpoint. Independent path integration gives 1.44549% and 2.98886%; the 240,000-step values differ by less than 0.0015 percentage point and converge toward them under refinement.

For transparency, the ESP metric is the root-mean-square difference in real $q/p_i$ at the digitized experimental values of $p\prime /p_i$, multiplied by 100. The computed curve is interpolated in pressure; this is not an orthogonal distance or a strain-matched residual. The common calculated locus can consequently agree with the same stress points while the two internal parameterizations differ. A physical flow comparison would need an independently defined strain mapping and strain-sensitive observations.

Figure 3 should be read in two stages. Panel (a) shows that the adopted surface, volume balance and hardening constrain both calculations to essentially the same curve in effective-stress space. Panel (b) then shows that the two calculations use different amounts of the internal coordinate to traverse that curve. These observations are complementary rather than contradictory: the first concerns the locus, and the second concerns its numerical parameterization. In particular, the approximately twofold difference between the marked $\chi$ values must not be translated into a twofold physical strain, time or stiffness without the missing kinematic and work-conjugate mapping.

5.2. Why Direct Overconsolidated Substitutions Need Further Equations

An overconsolidated initial state lies inside the current outer surface. The inherited bounding-surface architecture therefore introduces an image point on that surface and a radial mapping ratio $R$. In the supplied implementation, $p_{\mathrm{i}\mathrm{m}\mathrm{g}}=P_0\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(\tilde{\eta }/M\right)}^{\Psi }/\Omega \right]$ and $R=\mathrm{m}\mathrm{i}\mathrm{n}\left(p\prime /p_{\mathrm{i}\mathrm{m}\mathrm{g}},1\right)$. The ratios $M_d=MR$ and $M_k=M{R}^{-n}$ enter different equations: $D_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}=M_d-\tilde{\eta }$ controls plastic volume change, while $M_k-\tilde{\eta }$ controls the loading modulus. The inherited exponent is $n=0.5+2/\mathrm{O}\mathrm{C}\mathrm{R}$. This separation is important because an interior state can require a different loading stiffness even when the flow changes sign.

We examine two explicit choices. F1 replaces $D_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}$ by Equation (3) in the plastic increments but retains the inherited loading modulus. F2 also replaces the dilatancy-dependent factor in that modulus. With $c_p={\lambda }_c-{\kappa }_c$ and the positive common flow scale $g_{\tilde{q}}=1/\left(M{p}_{\mathrm{i}\mathrm{m}\mathrm{g}}\right)$, they are

$$K_p^{\mathrm{F}1}={\displaystyle \frac{\Omega }{c_p}}\left(M_k-\tilde{\eta }\right)g_{\tilde{q}}, K_p^{\mathrm{F}2}={\displaystyle \frac{\Omega }{c_p}}{\mathcal{D}}_{t_{*}}^{\alpha }g\hspace{0.17em}{g}_{\tilde{q}}, c_p={\lambda }_c-{\kappa }_c.$$

(15)

The candidate in these tests is evaluated at the actual normalized pressure $p\prime /P_0$, with the stated cutoff. The surface normal and common scale are evaluated at the image point. This choice is part of the experiment on equations, not a general fractional bounding-surface prescription. F1 and F2 are neither new validated models nor implementations of the complete architectures in [7,17]. Their purpose is to determine which inherited relations survive each direct substitution.

5.2.1. Analytical Sign Boundary for F2

Equation (6) is positive on one side of $t_{*}$ and negative on the other. Since the kernel in Equation (3) is positive, the integral has the same respective signs for every $0<\alpha <1$ and $\Omega >0$. At the isotropic overconsolidated initial state, ${\bar{p}}_0=1/\mathrm{O}\mathrm{C}\mathrm{R}$ and $t_0=\mathrm{l}\mathrm{n}\left(\mathrm{O}\mathrm{C}\mathrm{R}\right)$. All other factors in $K_p^{\mathrm{F}2}$ are positive, hence

$$\mathrm{sgn}{K}_p^{\mathrm{F}2}=\mathrm{sgn}(1/\Psi -\mathrm{lnOCR}),{\mathrm{OCR}}_c=\exp (1/\Psi ).$$

(16)

For $\Psi =1.20$, the boundary is ${\mathrm{O}\mathrm{C}\mathrm{R}}_c=2.301$. In F2, a candidate dilative response at a larger OCR is therefore tied directly to a negative loading modulus. Dilative plastic volume change is not itself the failure: the failure is the loss of the positive-modulus loading structure that this substitution is intended to retain. Figure 4a shows that changing $\alpha$ changes the magnitude but cannot move the boundary. The numerical illustration spans $1.001\le \mathrm{O}\mathrm{C}\mathrm{R}\le 20$; the exactly isotropic NC evaluation uses the stated cutoff.

5.2.2. Transformed Plastic Work and Consistency Conditions for F1

For the declared transformed-stress pairing, the plastic coordinate increments are $\mathrm{d}{\epsilon }_v^p=\mathrm{d}\Lambda \left({\mathcal{D}}_{t_{*}}^{\alpha }g\right)g_{\tilde{q}}$ and $\mathrm{d}{\chi }^p=\mathrm{d}\Lambda g_{\tilde{q}}$. During active plastic loading with $p^{\prime }\text{>}0$, $\mathrm{d}\Lambda \text{>}0$ and $g_{\tilde{q}}\text{>}0$, Equation (5) reduces to

$$\mathrm{d}{W}^p=d\Lambda \hspace{0.17em}{g}_{\tilde{q}}\hspace{0.17em}p\prime \left({\mathcal{D}}_{t_{*}}^{\alpha }g+\tilde{\eta }\right)\ge 0 \iff {\mathcal{D}}_{t_{*}}^{\alpha }g+\tilde{\eta }\ge 0.$$

(17)

For positive $p\prime$, $g_{\tilde{q}}$ and active plastic multiplier, the bracket is the entire sign condition for the declared work pairing. A negative candidate $D$ makes the volumetric contribution negative; sufficient positive shear work can compensate only when $\tilde{\eta }\ge -D$. Figure 4b shows that balance. When the plastic multiplier is zero, the assigned plastic work is zero irrespective of the bracket. This interpretation explains the computed sign without equating it to a general thermodynamic dissipation theorem.

The other F1 defect concerns the connection between volume change and hardening. At normal consolidation, consistency requires the stress increment and the surface-size increment to satisfy the differential of Equation (1). Because $P_0$ changes through the adopted plastic volumetric increment, its contribution to the loading modulus must use the same $D$ as the flow. Retaining the base factor after changing $D$ leaves

$$\Delta K_p=K_p^{\mathrm{F}1}-K_p^{\mathrm{consistent}}={\displaystyle \frac{\Omega }{c_p}}\left[\left(M_k-\tilde{\eta }\right)-{\mathcal{D}}_{t_{*}}^{\alpha }g\right]g_{\tilde{q}}.$$

(18)

The mismatch is generally non-zero, so the computed plastic multiplier does not enforce the surface consistency equation for the substituted flow. This is a structural difference between equations, rather than an increment-size error. Equation (18) is a necessary check at normal consolidation; for interior states, a complete overconsolidated derivation must also include the image-mapping structure. Appendix B makes explicit that no surface projection is applied to hide the residual.

The two panels diagnose different couplings. In panel (a), a negative $D$ represents dilation under the declared sign convention; dilation itself is physically plausible for overconsolidated clay. The F2 problem is narrower: that sign is transferred directly into a loading modulus that the inherited bounding-surface architecture intends to remain positive. Panel (b) asks instead whether the positive deviatoric contribution is sufficient to offset the negative volumetric contribution in the declared pairing. Neither curve is a yield surface, a failure envelope or an experimentally identified OCR boundary. They are sign boundaries for the two precisely stated trial substitutions.

Table 5 combines the analytical F2 initial-state sign test with integrated F1 diagnostics. No complete F2 overconsolidated trajectory is claimed. For F1, the maximum residual tests whether the stress has crossed outside the outer surface, while the minimum work rate tests the separately declared sign condition.

At OCR 6 and 10, F2 starts with negative $K_p$, in agreement with Equation (16). This violates the inherited positive-modulus requirement, but does not prove that every possible solver with a negative plastic modulus is singular. F1 remains numerically integrable over the reported runs. At OCR 1 and 1.2, however, its maximum $r_F$ is about 0.28: the computed stress has moved outside the surface whose hardening was prescribed. At OCR 6 and 10 the maximum residual remains negative, but the assigned plastic-work rate becomes negative. The two diagnostics therefore detect different problems. A dimensionless residual of 0.28 is not a 28% stress error.

Halving and doubling the increment count preserves these signs. The largest change in maximum residual is $8.9\times {10}^{-6}$, and the largest change in minimum work rate is below $1.6\times {10}^{-9}$ in the reported normalized units. The persistence under refinement supports the analytical diagnosis of the specified substitutions. It does not reject other fractional bounding-surface architectures, such as those in [7,17], which have their own gradients, surfaces and hardening laws. Earlier base-model fits cannot supply the missing validation for the different equations tested here.

6. Discussion

6.1. Connecting Each Mathematical Result to a Physical Question

The terminal calculation concerns the state at which the selected volumetric-flow mechanism ceases. The coordinate-composition calculation concerns whether the same scalar operation survives a change in stress representation. The incremental tests concern whether flow, pressure change and hardening still obey the equations declared for the embedding. These questions require different evidence. In particular, a mathematically unchanged integral cannot establish that its zero is the intended material terminal, and a coincident zero cannot establish consistent evolution toward it.

The extension/compression comparison makes the distinction concrete. For a transformed ratio $\tilde{\eta }$, the axisymmetric SMP inverse gives a real TE/TC ratio $3/\left(3+\tilde{\eta }\right)$. At $\tilde{\eta }=M$ this becomes $M_e/M$. At the incompatible operator zero, it must instead be evaluated at $\tilde{\eta }=M{\left(\Omega /\Psi \right)}^{1/\Psi }$. For kaolin, these are 0.7407 and 0.7825. The difference follows from evaluating the same stress transformation at different terminals, not from a failure of the transformation. Figure 5 shows how the separation depends on the two independent surface parameters.

The direction of the mismatch is systematic. When $\Omega /\Psi <1$, the operator terminal lies below the base terminal in transformed stress ratio, while its normalized surface pressure is higher; when $\Omega /\Psi >1$, both directions reverse. Only $\Omega /\Psi =1$ makes the two ratios equal to unity. The colors display these deterministic consequences of the two selected surface parameters. They do not represent confidence limits, data density or prediction error, and the pressure ratio in panel (b) remains distinct from the evolving normally consolidated mean-stress ratio reported in Table 3.

6.2. What the Data can and Cannot Identify

Under Equations (13)–(14), the NC stress locus contains no independent flow-direction information. Its RMS discrepancy cannot therefore identify $\alpha$. The numerical loading coordinate can distinguish the two embeddings, but an experimental comparison would first require the mapping between that coordinate and a measured shear strain, together with compatible elasticity. Figure 3b is useful as a controlled calculation of how the declared equations progress; it is not evidence that the clay needs twice as much physical strain under fractional flow.

The roles of the inherited datasets are correspondingly unequal. Fujinomori supplies stress-path context for the NC identity. Kaolin supplies coefficients for a terminal comparison and substitution diagnostics. Eastern Osaka supplies only an algebraic parameter illustration. No new fit or independent identification of fractional order is performed. The documented Fujinomori normalization in [12] should not be extended automatically to the other coefficient pairs, and preserving digitized files does not establish their original experimental accuracy.

6.3. Consequences for Further Constitutive Development

For the terminal issue, imposing $\Omega =\Psi$ would make the two selected zeros coincide, at the cost of removing one degree of freedom in the surface geometry. Retaining independent parameters would instead require a different operator terminal or a revised material critical-state specification. Either choice would require a new consistency derivation; the current results do not demonstrate that simply moving the integration limit preserves the other properties of Equation (3).

For physical validation, a complete stress–strain relation is needed before laboratory shear or volumetric measurements can distinguish the proposed flow. A useful experimental comparison would then include observables sensitive to contraction and dilation, with the loading path and terminal behavior documented. Overconsolidated and intermediate-Lode responses require their own loading, direction and hardening definitions. The present equations do not provide those missing parts merely by fitting additional parameters.

The immediate result is a set of constraints on model construction. The shared NC locus explains why further stress-path fitting alone cannot settle the flow question. The incompatible zero changes the candidate endpoint. The two direct substitutions reveal where hardening and the chosen work condition cease to match the adopted flow. These findings narrow the next task to completing and testing the constitutive equations.

6.4. Relation to Established Clay-Plasticity Architectures

Critical-state formulations assign different roles to elasticity, the loading surface, plastic flow and hardening [1,2,33]. Bounding-surface plasticity additionally defines an image rule and a modulus for states inside the outer surface [22,34,35]. Natural and anisotropic clay models may include fabric evolution, rotational hardening, non-coaxiality and loss of structure [23,36,37]. Calibration is conditional on that complete architecture and the selected observables [24]. Our F1/F2 results illustrate why a fractional replacement of one ingredient must be accompanied by a consistent definition of the others.

A physical work statement also requires care. Transformed stresses are a mathematical representation of the same material state, but an arbitrary stress–coordinate product is not automatically physical work. A successor formulation must specify a conjugate strain mapping and stored energy before it can use a dissipation inequality as a physical admissibility test. Here, the limited pairing in Equation (5) is stated explicitly so that its positive and negative values can be interpreted without a stronger claim.

7. Limitations

The analysis concerns a scalar phenomenological construction on an isotropic, rate-independent meridian. It provides no micromechanical derivation of the fractional kernel, time-memory law, complete tensorial strain mapping, hydrostatic regularization, consistent finite-element tangent or boundary-value prediction. Eastern Osaka is natural and structured, so its parameters enter only the conditional algebraic terminal comparison. F1 and F2 are specified diagnostic substitutions and do not exhaust fractional plasticity.

The available records are digitized from published figures, rather than raw laboratory files, and only the Fujinomori stress path is used as experimental context for a numerical curve in this study. Laboratory strain does not validate $\chi$. Fractional order is inherited rather than fitted, and the near-isotropic cutoff is part of the numerical definition. The conclusions apply to the stated coefficients, surface, hardening and work pairing; numerical reproduction cannot remove these physical and evidential limits.

8. Conclusions

The analysis was organized as a physical sequence: identify what the operator is intended to represent, determine whether its zero preserves the material terminal, establish whether the scalar operation survives the stress transformation, and test whether the resulting increments remain consistent with hardening and the declared work pairing. The results show that these requirements are logically separate.

(1) The terminal-directed integral fixes a zero of the candidate volumetric flow at the surface stationary point. That zero coincides with the non-associated base-model terminal if and only if $\Omega =\Psi$. For kaolin $\left(\Psi =1.20,\Omega =0.91\right)$, the alternative zero occurs at $\tilde{\eta }/M=0.794$ and gives a candidate NC undrained mean stress 18.6% above the base terminal. This is a change in the selected model endpoint, not a measured error.

(2) The scalar operator is invariant under a stress transformation that preserves its integration coordinate, terminal and kernel. SMP satisfies this sufficient condition for the present pressure-based coordinate. A generic nonlinear map gives 19.704% discrepancy for the ordinary local-chain substitute at $\alpha =0.641$. Neither result supplies the physical strain mapping needed by a complete constitutive model.

(3) The NC undrained stress locus follows from the surface, elastic volume balance and hardening without selecting a flow direction. Accordingly, the Fujinomori stress-path discrepancies are 4.312% and 4.313% for the two candidates, while the first 99% terminal crossings occur at about 1.45% and 2.99% of the numerical coordinate. The latter values demonstrate different calculated progression and are not laboratory-strain predictions.

(4) The direct OC substitutions expose two different inconsistencies. F2 couples the loading modulus to a sign-changing candidate, reversing its initial sign for $\mathrm{O}\mathrm{C}\mathrm{R}>\mathrm{e}\mathrm{x}\mathrm{p}\left(1/\Psi \right)$; the kaolin threshold is 2.301. F1 retains an incompatible hardening contribution and exhibits exterior surface drift at low OCR or negative work under the declared pairing at high OCR. Refinement preserves the reported signs.

(5) A fractional flow construction must connect its selected zero, stress representation and incremental hardening to the intended physical response. The present study establishes constraints on that connection. Experimental validation of a completed model requires strain and volumetric observables, a defined stress–strain mapping and a justified physical dissipation statement.

The practical implication is that additional stress-path fitting is not the next constitutive step. A successor model should first select and justify its material terminal, derive the loading modulus and hardening law consistently with the selected flow measure, and formulate a strain mapping conjugate to the transformed stress. It can then be tested against measured shear and volumetric strains across OCR and Lode angle, with an explicit stored-energy statement where a physical dissipation inequality is used. Only after those links are specified can the fractional order be identified from physical observations rather than retained as an inherited numerical setting.