Long-form technical specification.

0A.1 Coverage

The documentation covers the current implemented workflow: CPT file parsing for GEF, Excel, and CSV sources, per-point CPT classification, layer detection, layer parameter assignment, stiffness derivation, experimental m-fitting, and the Stage 6 engineering applications.

Where the application contains an engineering simplification, the simplification is stated explicitly. The purpose is not to replace engineering judgement, but to make the implemented mathematics auditable and readable.

• Classification and boundary logic are documented separately from parameter assignment.
• Stage 6 uses the interpreted CPT state produced by the earlier stages rather than reclassifying the profile independently.
• The notation below follows geotechnical convention with compression positive and depth measured downward.

0A.2 Layer quantities carried into engineering use

Per interpreted layer, the application carries the quantities needed by the engineering modules: top and bottom depth, γ and γ_sat, φ′, c′, c_u, E_oed,ref, E_oed,i, E_50,ref, E_ur,ref, m, K_0,nc, ν_ur, k_h, and k_v.

Geometry, load combination, hydraulic scenario, and optional Stage 5 tuning are then applied on top of those layer values.

1.1 Supported file structures and column mapping

GEF files declare their physical quantities through #COLUMNINFO lines. The application reads the quantity identifier and maps it to the relevant column index; column order is therefore never assumed.

Excel workbooks are read from a Data sheet and, when present, a Header sheet. The Data sheet may contain depth, q_c, f_s, and R_f; the Header sheet can provide project, test, location, date, water level, surface level, coordinates, and the net area ratio.

CSV files are the reduced-input route. They must contain headers named depth and qc. Optional headers are fs and rf. Comma, semicolon, and tab delimiters are detected automatically.

A recommended minimal CSV header is depth [m], qc [MPa], fs [MPa]. If decimal commas are used in the values, a semicolon-delimited CSV is recommended.

• GEF quantity 1: penetration length.
• GEF quantity 2: q_c.
• GEF quantity 3: f_s.
• GEF quantity 4: R_f.
• GEF quantity 6: u₂.
• GEF quantity 11: corrected depth, used before quantity 1 when both are present.
• Excel/CSV column headers may include units, for example qc [MPa], fs kPa, or Friction ratio (Rf) in %.

1.2 Unit conversion and row filtering

The parser converts q_c and f_s to MPa based on the declared GEF unit string or the Excel/CSV column header. Headers containing MPa are used directly; kPa is divided by 1000; Pa is divided by 1 000 000.

If the unit declaration is absent or ambiguous, a heuristic fallback is used so the file remains workable: q_c values above 100 are treated as kPa, otherwise MPa; f_s values above 1000 are treated as Pa, values above 10 as kPa, otherwise MPa. Explicit unit labels are therefore strongly recommended.

Rows are removed when they clearly do not represent meaningful CPT measurements, namely negative depths, values before cone engagement, or all-zero terminal rows.

• Depth is expected in metres below surface.
• R_f is optional; when absent or invalid, the app computes it from f_s and q_c.
• Rows with z < 0 are discarded.
• Rows with q_c < 0.02 MPa are treated as cone-not-engaged and discarded.
• All-zero rows appended by logging software are discarded.

1.3 Water table, elevation, and header metadata

The phreatic level is first taken from the GEF measurement variables or the Excel Header field Waterniveau when available; otherwise a default depth below surface is assigned. CSV files do not carry this metadata, so the engineer should review or override the default after import.

Surface elevation is read from the GEF ZID header, the Excel Header field Grondniveau, or entered manually. When present, all depth values can be expressed relative to TAW as well as depth below surface.

2.1 Robertson (1990) — normalized SBT / I_c

The Robertson route uses a normalized q_t–F_r framework. Because final layer unit weights are not yet available at Stage 2, the app uses a preliminary stress estimate with fixed unsaturated and saturated unit weights.

If pore pressure u₂ and net area ratio a are available, the cone resistance is corrected to q_t; otherwise the measured q_c is used directly.

• The current app groups the normalized result into Gravel, Sand, Silty sand, Sandy clay, Clay, and Peat / organic.
• Sensitive fine-grained soil is not inferred separately from I_c alone.

2.2 Robertson (2016) — iterative normalized SBT / Q_tn

The Robertson 2016 route keeps the same broad SBT / I_c family mapping as the earlier normalized Robertson workflow, but replaces the fixed-stress-exponent cone normalization with the iterative Q_tn formulation. The stress exponent n varies with the inferred soil behaviour and is solved iteratively.

This route may be preferred when the input is a CPTu, but the implementation does not require measured u₂. If u₂ is absent, the app uses q_t = q_c in the same way as the Robertson 1990 fallback.

• The implementation reuses the existing broad app families: Gravel, Sand, Silty sand, Sandy clay, Clay, and Peat / organic.
• The displayed method metric is Q_tn; the I_c boundaries remain the same as the Robertson 1990 route.

2.3 CUR 3 layers — broad q_c–R_f zoning

The CUR 3 layers route is implemented as a direct q_c–R_f zoning rule based on the published broad chart with four material fields: Sand, Silt, Clay, and Peat. It is a boundary-generation method, not a detailed parameter catalogue.

The implementation follows the chart as nested zones checked from the most specific region to the broadest envelope: sand first, then silt, then clay, with peat as the remaining outside region.

To keep the downstream parameter workflow stable, the chart field “Silt” is carried internally as the app’s intermediate family and tagged with the subtype marker CUR3 silt. The chart logic itself remains four-zoned.

Show published CUR 3 chart

• This route is direct and non-normalized: no stress correction is applied to q_c before classification.
• The current app uses this route for broad layer boundary logic; detailed geotechnical parameter assignment is still handled later by the Stage 3 parameter method.
• Internal mapping note: the chart field “Silt” is stored as the app’s intermediate type family so that later parameter and compatibility tables remain usable.

2.4 NEN 6740 — stress-dependent material classification

NEN 6740 uses a stress-corrected cone resistance together with friction ratio. The source method is graphical: a semilog chart with fourteen material areas rather than a closed algebraic decision tree.

The app therefore separates the method into two parts. First it evaluates the published stress correction documented in the Deltares D-SHEET Piling manual. Then it applies a transparent representative-area rule that selects the nearest material area from a digitized fourteen-material set.

Show published NEN 6740 chart

• The stress-correction equation with exponent 0.67 is part of the documented Deltares implementation of the published NEN route; the one-dimensional score projection is the app’s deterministic way of reproducing the graphical chart.
• The score coefficient magnitude 0.34 is not a normative NEN coefficient. It is the regression-fit semilog projection magnitude through the fourteen stored representative centres used by the app.
• The representative materials are: gravel, slightly silty, moderate; sand, clean, stiff; sand, slightly silty, moderate; sand, very silty, loose; loam, very sandy, stiff; loam, slightly sandy, weak; clay, very sandy, stiff; clay, slightly sandy, moderate; clay, clean, stiff; clay, clean, weak; clay, organic, moderate; clay, organic, weak; peat, moderately preloaded, moderate; peat, not preloaded, weak.
• The weakest transition remains the area-5/area-6 pair (loam, slightly sandy, weak versus clay, very sandy, stiff). Near ties are resolved by the fixed chart order of the stored representative set.
• A repo verification script re-projects all fourteen centres at several effective stresses to ensure the stored centres and the implemented score formula stay in sync.
• The selected NEN material is then mapped into the app’s internal type families so that later layer grouping and parameter workflows remain consistent.

2.5 NEN Tabel 3 — direct subtype classification

The Table 3 route uses raw q_c and R_f and returns a detailed subtype together with the characteristic parameter row used later in the workflow. The direct source documented here is the implemented NEN Tabel 3 catalogue itself, not a separate Eurocode classification chart.

Because the table contains overlapping q_c–R_f envelopes, the app evaluates the rows deterministically in table order rather than assuming the source is mutually exclusive. The implemented family order is grind, zand, leem, klei, then veen. Within each family, subrows are checked top-to-bottom, with q_c lower bounds inclusive and upper bounds exclusive. For R_f, the open bands R_f < 1 and R_f > 6 remain strict, while the bounded intervals are treated as closed.

For this method, the raw boundary logic follows the detailed subtype result rather than only the broad family. This means subtype changes such as density or consistency transitions inside one family can create provisional layer boundaries before smart merge and minimum-thickness correction simplify them.

• If no table row matches, the app keeps a deterministic fallback so the workflow can continue, but that fallback is not itself part of Table 3.
• The app UI may still show a legacy “EC7” label because the later parameter workflow is Eurocode-aligned. The direct q_c–R_f table documented in this section is the implemented NEN Tabel 3 source.

2.6 Classification versus parameter assignment

Stage 2 is boundary logic only. The classification method decides which soil type each CPT reading belongs to and therefore where soil-type changes occur with depth.

Stage 3 parameter assignment is independent of the chosen classification route. A layer identified by Robertson, CUR 3 layers, or NEN 6740 may still be assigned a Eurocode Table 3 subtype and its associated γ, γ_sat, φ′, c′, and c_u values.

3.1 Raw segmentation and boundary placement

The first pass creates provisional raw segments directly from the pointwise Stage 2 classification sequence. For Robertson, CUR 3 layers, and NEN 6740, the raw segment key is the broad internal soil family. For NEN Tabel 3, the raw segment key is the pair {type, subtype}, so detailed subtype changes can create provisional layer boundaries before later simplification.

Boundary depths are then constructed from the retained CPT rows. The first raw segment starts at ground level. Other raw segments start at the midpoint between the last CPT row of the previous segment and the first CPT row of the new segment, unless a downward merge has already handed them an inherited top boundary. If that straddling sample pair is unavailable or degenerate, the app falls back to the legacy 20 mm offset from the first row of the new segment. Segment thickness is always evaluated on these active layer boundaries rather than on the bare difference between first and last CPT reading.

• This means the NEN Tabel 3 route is intrinsically finer than the other Stage 2 routes, because subtype transitions inside one soil family are preserved in the raw segmentation.
• NEN 6740 currently uses the internal family only for raw layer generation; its 14 detailed material labels remain available for interpretation but do not by themselves create raw boundaries.
• The midpoint rule is self-calibrating to the actual CPT logging step and remains bounded by construction between the two samples that straddle the raw boundary.
• For uniform 40 mm sampling, the midpoint rule collapses exactly to the earlier 20 mm offset.

3.1A Smart merge correction

The smart-merge correction is an in-house MADEP post-processing heuristic. It is not a published external classification standard or code method.

The app also provides a smart-merge mode. This mode is applied only after the original raw layering has already been created from the unmodified point-by-point classification sequence.

In smart mode, the algorithm first removes highly compatible existing boundaries when the continuity score exceeds a sensitivity-controlled threshold. Only after that reduction step does it enforce the minimum thickness by merging every remaining sub-minimum layer with its most similar neighbor.

The resulting smart-merged configuration becomes the active layer model. Layer boundaries, per-layer averages, subtype suggestion, and downstream parameter assignment are recalculated from the final merged geometry.

• The original boundary detection step is never redefined by smart merge; smart merge operates only after the base layering algorithm has already run.
• In smart mode, similarity-based reduction is applied before the final minimum-thickness enforcement. The minimum thickness therefore acts as the last hard constraint, not as the first distortion of the layering.
• Very thin sections are intentionally down-weighted in the resistance to merging: sliver layers receive a small merge bonus, and sharp-transition penalties are attenuated when the layer thickness is small relative to a fixed sliver reference. This early similarity reduction is therefore independent of the chosen minimum thickness.
• Low sensitivity keeps smart merge closer to the conservative baseline behavior; high sensitivity allows smaller continuity advantages to influence the merge direction and makes it easier to remove highly compatible boundaries. At the extreme right-hand end of the range, the smart pass becomes willing to absorb almost all remaining compatible boundaries before the minimum-thickness rule is applied.
• The final layer count is not strictly monotonic with the sensitivity value. Because minimum-thickness enforcement runs after the smart-merge pass, a higher sensitivity can change the topology of the intermediate layering in a way that later produces either fewer or more final layers than a lower sensitivity case.
• If the two candidate directions are nearly equal, the app resolves the tie by merging into the thicker neighboring layer; if the thicknesses are equal, upward merge is chosen.
• The penalty term reduces the score for large q_c jumps, large R_f jumps, peat/non-peat transitions, gravel/non-gravel transitions, and thin layers that act as critical markers such as very weak peat, very coarse gravel, very high R_f, or extreme q_c.

3.1B Minimum-thickness enforcement

The minimum layer thickness remains the final hard constraint on the smart-merge chain. In baseline mode, it is enforced directly by repeated upward merging. In smart mode, it is enforced only after the similarity-based boundary reduction described above.

In smart mode, each remaining sub-minimum layer is merged with the more similar neighboring layer according to the directional score.

• This means the minimum-thickness setting does not control the early smart-similarity reduction itself; it only controls the last mandatory cleanup step.
• As a result, low t_min preserves more of the similarity-reduced layering, while high t_min forces stronger final consolidation of the remaining thin layers.

3.1C Intended engineer workflow

The intended use of the layering controls is sequential. First select the Stage 2 classification method that best represents the CPT behaviour for the project. That method defines the provisional pointwise geology and therefore the raw segment sequence.

Then use smart merge to remove boundaries that are judged highly compatible in the context of the surrounding profile. Only after that should the minimum thickness be used as the final hard cleanup rule for layers that remain too thin to retain.

The minimum-thickness setting is therefore not intended to replace the geological interpretation performed by the selected classification method or by the smart-merge pass. It is the final practical consolidation rule applied after those interpretation steps.

• Choose the classification method first; it defines the raw boundaries that everything else works from.
• Use smart merge to smooth obviously over-segmented or highly compatible boundaries, not to force an arbitrary target layer count.
• Use minimum thickness last to eliminate remaining thin layers that are still impractical after the similarity-based cleanup.
• After every change, review the resulting boundaries against the q_c and R_f trends rather than judging only by the final number of layers.

3.2 Final layer statistics and boundary continuity

After the final merged geometry is fixed, the app recomputes each layer summary from the merged CPT rows. Mean q_c, mean f_s, and mean R_f are taken over the valid CPT rows inside the final segment. The representative subtype string is the modal subtype within that final row set.

The final layer table is then rebuilt with continuous boundaries: each layer top is taken from the previous final layer bottom, and each layer bottom is forced not to lie above its own top. These recomputed averages are the Stage 4 input for stiffness and conductivity derivation.

• Rows with q_c ≤ 0.02 MPa are excluded from the averages when possible; if that would leave no rows, the full segment row set is used as fallback.
• So every time the merging changes, the boundaries and the engineering averages are recalculated consistently from the new merged geometry.

3.3 Subtype suggestion and parameter assignment

The final merged layer first inherits provisional parameters from the merged rows. For the NEN Tabel 3 classification route, the app can average the pointwise row parameters available from the table-matched CPT rows. For the other classification routes, it falls back initially to the generic family defaults.

After that, the app auto-suggests the best Eurocode Table 3 subtype from the catalogue using the final layer average q_c, average R_f, and the compatibility matrix of the CPT family. If the active Stage 3 parameter method is the Eurocode Table 3 route, the suggested subtype also supplies the layer parameters. All assigned values remain editable, and manual overrides always take precedence over automatic assignment.

• The suggested Eurocode subtype is selected from the full catalogue by compatibility first and q_c/R_f fit second.
• The final layer warnings below the table are then derived from the chosen subtype versus the CPT-family compatibility matrix, so the engineer can still see when a selected subtype lies only in an adjacent transition family.

3.4 PLAXIS simulated CPT export

Stage 3 can export the active interpreted layer model as a measurement-style CPT text file for PLAXIS. The goal is to reproduce the final interpreted layering exactly, not to retain the original pointwise cone variability.

The export reuses the original CPT depth rows and coordinates. For each original depth, the app finds the active final layer and writes that layer's representative values back as a synthetic CPT row. The result is therefore piecewise constant per interpreted layer, but sampled at the original CPT spacing.

• The header contains X[m], Y[m], Z[m], followed by a D[m] Q[MPa] F[MPa] x table.
• One exported row is written for every original CPT depth row, so the synthetic file remains aligned with the measured depth sampling.
• The final x column is currently written as 0, so the exported file behaves as a PLAXIS-compatible CPT with depth, cone resistance, and sleeve friction.

4.1 Effective stress at layer midpoint

Both model-parameter routes use the effective vertical stress at layer midpoint. The phreatic level therefore directly affects the reference stiffness correction.

Above the water table, γ is used; below it, γ_sat is used. Pore pressure is hydrostatic.

4.2 q_c-to-E_oed,i correlation and α methods

The first stiffness step converts the representative cone resistance into an oedometric modulus through one of two α correlations selected in the Stage 4 UI: Method A — Sanglerat literature (fixed α per behavioural soil type) or Method B — SB260-21-6.4.10 (q_c-graded family rules). Method A and Method B are alternative routes that produce different α — and hence different E_oed,i — for the same layer.

Method A applies fixed Sanglerat-style α defaults keyed to the broad layer type (Peat, Soft clay, Clay, Sandy clay, Silty sand, Sand, Gravel). All rows in the layer share one α regardless of q_c.

Method B applies the SB260-21-6.4.10 Tabel 21-6-5 rules, structured around three families — cohesive (klei, leem, veen), transition / overgangsgronden (the (zh) and (lh) subtype variants), and granular / zandgronden (zand, grind, grind (kh)). The family is selected from the EC7 subtype string first, falling back to the broad type when no subtype is set, and α (or E_s) varies with q_c within the family.

The full implemented α mapping is shown in the expandable reference table below so the reader can audit exactly which family, q_c band, or modulus formula is applied.

• Method A produces the higher α values for granular soils (Sand = 13, Gravel = 15) characteristic of Sanglerat literature. Method B uses the SB260 zandgronden rule E_s = 4q_c / 2q_c+20 / 120, which gives much smaller effective α (e.g. α = 4 for q_c ≤ 10 MPa). The two methods are not interchangeable — Method A and Method B will give different E_oed,i for the same layer.
• Method A keys off the broad layer type (l.type) only. Method B keys off the EC7 subtype string first and falls back to type when no subtype is available.
• For peat, water content w is not available in the app, so Method B uses the SB260 default α = 1.5 (the conservative branch of Tabel 21-6-5 for veen).
• The (zh) / (lh) suffix on klei / leem / zand routes to the transition family. The (kh) suffix on grind routes to the granular family.

4.3 Reference-stress correction (shared by Methods A and B)

Both Method A and Method B apply the same Hardening Soil reference-stress correction to convert the in-situ E_oed,i from §4.2 into the reference-stress quantity E_oed,ref at p_ref = 100 kPa. The cohesion-corrected form is used so the correction remains well-behaved in cohesive layers where σ′_v0 can be small relative to c′cotφ′.

Methods A and B differ only in how E_50,ref and E_ur,ref are derived from this shared E_oed,ref; the depth correction itself is identical in both routes.

The CUR 2003-7 stress-exponent default is binary: m = 0.5 for granular soils (Sand, Silty sand, Gravel) and m = 1.0 for cohesive soils (Sandy clay / leem, Clay, Soft clay, Peat / organic). Stage 5 m-fitting can override the m default per layer.

4.4 Method A — CUR 2003-7 stiffness ratios

Method A takes the shared E_oed,ref from §4.3 and applies the CUR 2003-7 family rule to derive E_50,ref. E_ur,ref is taken as three times E_50,ref. The family rule treats klei and leem together, so Sandy clay (leem) is in the cohesive set with the 1.25 factor; for granular soils E_50,ref is set equal to E_oed,ref.

4.5 Method B — E_50,ref = E_oed,ref

Method B takes the shared E_oed,ref from §4.3 and sets E_50,ref equal to E_oed,ref for all soils. E_ur,ref remains three times the selected E_50,ref.

This gives a single consistent reference stiffness and is sometimes preferred in practice when the engineer wants to avoid the cohesive-soil E₅₀/E_oed split.

4.6 Hydraulic conductivity basis

The app uses indicative hydraulic conductivity values tied to Belgian and USDA-style texture classes, with OVAM and De Smedt as the principal reference sources. The representative value is treated as a geometric-mean estimate within the adopted class range rather than a deterministic measurement.

Anisotropy is then introduced through a k_h/k_v ratio: isotropic for coarse granular soils, higher horizontal-than-vertical conductivity for fine-grained and layered soils.

For the PLAXIS material-command export, the internally stored conductivities are converted from m/s to m/day before they are written to the command file.

• Indicative anisotropy in the current logic is 1 for sand and gravel, and 3 for clay, loam, silty soils, and peat.
• In-situ measurement takes priority over indicative table values.

4.7 PLAXIS material-command export

Stage 4 includes a dedicated PLAXIS material export in addition to the CSV layer table. The implemented export writes a text file containing supported soilmat commands, not a directly generated native .matXdb database.

That choice is deliberate: in the modern PLAXIS workflow the app reliably creates project materials through command-line material creation, after which the engineer can copy those project materials into a reusable global material database from inside PLAXIS if desired.

For every interpreted layer, the export writes two material definitions: one Mohr-Coulomb material and one Hardening Soil material.

• Mohr-Coulomb export writes Identification, SoilModel = 2, DrainageType, gammaUnsat, gammaSat, ERef, nu, cRef, phi, psi, PermHorizontalPrimary, and PermVertical.
• Hardening Soil export writes Identification, SoilModel = 3, DrainageType, gammaUnsat, gammaSat, E50Ref, EOedRef, EURRef, PowerM, pRef = 100, cRef, phi, psi, PermHorizontalPrimary, and PermVertical.
• For the MC material, ERef follows the current-stress loading stiffness E50,i from the selected Stage 4 stiffness route rather than the Hardening Soil reference-stress quantity E50Ref.
• Method A gives E50,i = 1.25 EOed,i for Clay, Soft clay, and Peat / organic, and E50,i = EOed,i for the other soils. Method B gives E50,i = EOed,i for all soils.
• Only clean sand and clean gravel are exported as Drained by default. Fines-bearing granular subtypes such as zand (lh) and grind (kh) are exported as Undrained A.
• The export intentionally omits RF, nuUR, and K0NC in the current tested PLAXIS workflow so PLAXIS can keep those automatic or read-only values under its own control.
• cu and psi_unsat are not written by the material-command export. Undrained A uses the effective stress strength parameters, and psi_unsat remains available in the app without yet being exported.

5.1 Log-linear basis of the fit

The oedometer law of the Hardening Soil model becomes linear in logarithmic form. This allows the app to fit a straight line to the pointwise CPT-derived E_oed,i values inside one layer.

The fitted slope is interpreted directly as m, while the fitted intercept gives E_oed,ref.

5.2 OLS solution, fit quality, and acceptance

The regression is performed on all valid CPT points inside the layer, using the current α method and the current water table. The resulting m is clamped to the engineering range used by the app, and the fit quality is reported in log space.

The current interface shows the default line, the fitted line, and the point cloud; the engineer may still adjust the previewed m before acceptance.

• The fit is experimental and remains an engineer preview.
• m is expected to remain between 0 and 1 in normal interpretation, even though the raw regression may suggest otherwise.
• Accepting the fit overrides m only; the reference stiffness is then recomputed from the accepted m.

6.1 Sign convention and in-situ stresses

Compression is taken as positive. Depth z is measured downward from ground level. Elevation values are interpreted in Belgian TAW convention where relevant.

The in-situ effective-stress profile is reconstructed from the active groundwater level and the interpreted unit weights.

• Above the phreatic level, the dry unit weight γ is used; below it, γ_sat is used. Pore pressure u(z) = γ_w · max(0, z − z_w) is subtracted from σ_v(z) to form σ′_v(z).

6.2 Hardening Soil reference-stress convention

Stage 6 uses the Stage 4/5 Hardening Soil stiffness convention. E_oed,ref is referenced at p_ref = 100 kPa and stress-dependent stiffness is recovered from the interpreted m-value and effective stress level.

For the current settlement and dewatering routes, the governing stress is the vertical effective stress on the evaluated centreline.

6.3 Integration and sublayering

The engineering calculations use stratification summation. The profile is divided into sublayers, stresses are evaluated at sublayer mid-depth, and settlement or stress changes are integrated over depth.

The application keeps interpreted layer boundaries and phreatic boundaries explicit in the discretisation so stiffness and stress changes are not smeared across those interfaces.

7.1 Current implemented resistance model

The application computes drained and undrained ultimate resistance separately and converts those results to q_d or q_allow depending on the selected safety route.

The underlying reference family is the Brinch Hansen / EC7 Annex D shallow-foundation format. The current app implements the vertical-load subset of that model: no horizontal load, no ground slope, and no base tilt. Eccentricity is only used to derive the effective dimensions for the shape-factor route.

• The app now keeps the EC7 Annex D rough-base N_γ = 2(N_q − 1)tanφ′ formulation consistently. The earlier Vesić variant 2(N_q + 1)tanφ′ is slightly larger; Annex D is now preferred here because Stage 6 is an EC7-oriented screening tool and the Annex D form is slightly more conservative and easier to audit against the code text.
• Shape factors now default to Brinch Hansen / Annex D values derived from the effective plan ratio r = B′/L′. Users can switch to a conservative mode with all shape factors fixed to 1.0.
• For circular footings screened through this rectangular interface, use B = L and keep e_B = e_L = 0 so r = 1.
• For φ′ = 0, the drained block collapses to the undrained Prandtl expression q_ult,u = q + 5.14c_us_cud_cu; the drained expression is therefore not evaluated in that limit.
• The bearing-depth envelope is sampled from the selected starting depth upward to the model bottom using a depth increment clamped between 0.10 m and 0.25 m.
• This remains a shallow-foundation screening model: full inclined/eccentric-load verification, sliding, ground-slope effects, base tilt, and the full three-case groundwater averaging rule for the N_γ term are still outside the current app.

7.2 Belgian DA1 handling in the current app

For Belgian EC7 practice, the application evaluates DA1/1 with M1 soil strengths and DA1/2 with reduced M2 soil strengths. It then forms the governing drained and undrained envelopes separately.

The current Stage 6 app is therefore a resistance-side Belgian screening tool rather than a full action-side ULS verification engine.

• DA1/1 uses unfactored M1 soil strengths.
• DA1/2 uses reduced M2 soil strengths.
• The governing drained and undrained results are formed independently.

8.1 Scope and radius of influence

The hydraulic part of the module estimates drawdown at the CPT for a single well, an equivalent-radius excavation, or a line dewatering trench. The engineering deliverable is not the pumping rate alone, but the stress change and settlement induced at the CPT.

The current app still uses Sichardt to screen the extent of influence. This remains explicitly a rule-of-thumb estimate rather than a rigorous groundwater-flow boundary.

• C = 3000 is the classic Kyrieleis & Sichardt coefficient used in sandy-soil rule-of-thumb practice.
• Louwyck et al. (2022) show why Sichardt should be treated as a screening estimate, not as a rigorous groundwater-flow boundary.

8.2 Steady-state inflow — transmissivity-based screening model

The model treats the pumped zone as a stack of horizontal layers, each with its own horizontal conductivity k_h. For any chosen saturated thickness, the saturated part of each layer contributes to the total transmissivity according to its own thickness and conductivity.

Instead of representing the whole profile by one fixed conductivity, the model first constructs the transmissivity of the saturated profile and then updates that transmissivity as the phreatic level moves through the layered ground.

For unconfined flow, this changing hydraulic capacity is handled through the cumulative transmissivity moment. In that way, the drawdown solution remains analytical while still reflecting the layered structure of the interpreted CPT profile.

• T is the total horizontal flow capacity of the currently saturated part of the interpreted profile.
• T(h) expresses how that flow capacity changes as the water table rises or falls relative to the aquifer base.
• M(h) is the cumulative transmissivity moment used to solve the radial unconfined flow relation.
• For a homogeneous aquifer, the formulation reduces exactly to the classical Dupuit expression.
• The CPT does not measure b_i, T(h), or M(h) directly. The app derives them from the interpreted layer model: each layer receives k_h in Stage 4.6, the aquifer base is set by the hydraulic screening setup, and b_i(h) is then the overlap of layer i with the currently saturated interval between that base and the phreatic level.
• Once the interpreted layers and their k_h values are known, T(h) follows by summing k_h,ib_i(h) over the intersected layers, and M(h) is the corresponding cumulative transmissivity moment of that piecewise profile.

8.3 Drawdown profile toward the CPT

For radial unconfined flow, the app solves the head profile from the transmissivity-moment equation. For homogeneous conditions, this collapses to the classical Dupuit h² law.

For trenches, the displayed profile toward the CPT remains a linear screening interpolation, while the flow estimate itself is based on transmissivity.

8.4 Effective stress and settlement response

Once the new phreatic level at the CPT is known, the app recomputes pore pressure and effective stress with two selectable total-stress assumptions: conservative σ_v fixed, or realistic γ_sat → γ between the old and new water levels.

Settlement is then computed with the same constrained-modulus philosophy used in the settlement app, evaluated at the mean stress state between before and after drawdown.

Optional time curve for fine-grained layers follows the same Terzaghi 1D consolidation route used in the settlement application, with the consolidation coefficient c_v = k_vE_oed / γ_w derived from the interpreted layer k_v and the stiffness used for settlement.

• The app reports both S_conservative and S_realistic to expose the sensitivity to the chosen total-stress assumption.
• The public engineering output is total settlement versus distance from the source, with the active CPT position marked on that curve.

9.1 Settlement method summary

The current implementation follows the classical constrained-modulus route: compute q_net, derive Δσ_v beneath the foundation, evaluate E_oed at the mean stress level, and integrate ΔS = (Δσ_v / E_oed)Δz.

The evaluated location is the strip centreline for strip geometry and the centre of footprint for rectangular, square, and slab geometry.

9.2 Vertical stress increase — exact forms

For strip footings, the application can use the exact Boussinesq centreline solution. For rectangular loaded areas, the centre stress is computed by four-quadrant superposition of the corrected Newmark/Fadum corner influence factor.

The implementation includes the V < V₁ branch correction so the formula remains valid in the shallow, wide-load regime.

9.3 Constrained-modulus integration

For each sublayer, the app forms the mean effective stress between the in-situ and loaded state, evaluates E_oed at that level, and computes vertical strain and settlement increment from the constrained modulus.

This is the current implemented route for settlement in both sands and clays.

9.4 Output form and truncation

The review documents several truncation rules, but the current app lets the engineer choose the practical truncation setting. The present default in the interface is CPT bottom.

The output is a single vertical settlement beneath the evaluation point, not a 2D settlement field or edge-settlement map.

• Selectable truncation settings: Δσ_v < 10%σ′_v,0; Δσ_v < 20%q_net; CPT bottom.
• Optional time curve for fine-grained layers follows the same Terzaghi 1D consolidation route used in the dewatering application.

10.1 Modulus of subgrade reaction from CPT stiffness

The current implementation derives k_s from CPT-linked stiffness using the Vesić route. The app does not read E_s directly from a separate soil table at Stage 6; it derives the stiffness from the interpreted CPT layer model.

For each numerical sublayer between D_f and D_f + z_influence, the app evaluates E_oed at the local effective stress using the current layer stiffness law. It then converts that sublayer stiffness to the selected E_s route, averages it over the influence depth, and only then computes k_s. The result is therefore a footing-dependent support parameter, not a pure soil constant.

10.2 Governing equations and characteristic length

The beam is solved on either a Winkler or Pasternak elastic foundation. The corresponding differential equations are implemented numerically along the strip length.

The characteristic length determines whether the strip behaves as short, intermediate, or long on elastic support.

10.3 Pasternak 1D implementation

The Pasternak extension currently implemented in the app is a 1D strip formulation. The shear-layer parameter is not measured directly; it is inferred from the averaged soil shear modulus and a chosen influence depth.

In the current app, G_s,avg is obtained from the same CPT-derived stiffness profile used for k_s. The app first builds E_s,avg over the selected influence zone, then converts that average stiffness into the average shear modulus, and only then forms G_p.

This is why the Pasternak route is labeled as a screening extension rather than a continuum-calibrated Belgian design model.

• With the default app route, E_s,i = E_oed,i and ν_s = 0, so G_s,avg = E_s,avg/2.
• Uniform full-length loading can still produce very low bending moment because the deflection field approaches nearly uniform settlement.
• Patch loads and point loads are more informative when the objective is bending-driven reinforcement screening.

11.1 Structural design route from M_Ed

Once the ULS moment is obtained from the strip-on-foundation solve, the app converts concrete and steel to design strengths and estimates the required reinforcement area per meter width.

The current implementation also applies the EC2 durability and cover route, so c_nom is not just a free guessed number.

• The result should be read as strip-based ULS reinforcement screening, not as a final 2D slab reinforcement layout.
• Exposure class, working life, and detailing assumptions feed the EC2 cover route used in the calculation.

12.1 Present solver scope and geometry model

The current seepage module solves a steady-state Darcy/Laplace problem on a constrained triangular mesh. It is intended for cross-section screening, flow visualisation, exit-gradient checks, and optional pore-pressure handoff to Bishop.

Unless custom polygons are enabled, the soil layout still starts from the active CPT column extended laterally across the section. If custom polygons are enabled, those polygons become the seepage regions on the shared canvas.

• Outer-boundary edges carry Dirichlet, Neumann, and seepage-face conditions: terrain, base, and the lateral sides.
• Interior drain polylines carry head-prescribed Dirichlet with three gating modes (always, when-saturated, head-cap). The head-cap mode is a one-way head cap solved as a primal–dual active-set semismooth Newton loop on the LCP h ≤ H_d, R ≤ 0, (h − H_d)·R = 0, run jointly with the outer seepage-face active set.
• Other interior region boundaries remain meshing constraints only; they are not seepage BCs.
• The current module is steady-state only: no transient storage, no consolidation time response, and no Richards-type unsaturated flow law.
• Retaining walls are thin low-permeability regions with per-wall anisotropic conductivity (k_⊥, k_∥) in the wall-local frame. The current route supports axis-aligned vertical walls; general-orientation walls and true zero-thickness duplicate-node cutoffs are deferred.
• Flow-prescribed line sinks (q = −Q δ(x − x_w)) and finite-resistance Cauchy drains (q = β(h − H_d)) are not yet exposed.

12.2 Free surface, boundary logic, and numerical safeguards

The shipped default is the iterative free-surface mode. It keeps one fixed mesh, classifies elements as wet or dry from the current pressure-head state, and updates the atmospheric seepage-face activation set on the outer boundary.

A second fixed phreatic mode remains available for benchmarking or sensitivity studies. In that mode the drawn phreatic line is imposed as an internal constraint and the solution is locked to that water-table geometry.

• A prescribed-head BC on a sloping or vertical edge applies only to the submerged part below the entered head elevation; the dry part above reverts to natural no-flow.
• Iterative convergence is accepted only when the seepage-face active set is stable, the boundary flow-rate error is within the target tolerance, and no disconnected wet islands remain.
• The current default flow-rate error target is 1%, with a 10 s runtime cap in iterative mode.
• If the wet/dry loop creates isolated wet pockets that are not connected to a prescribed-head boundary, the solver applies a connectivity fallback and dries those pockets out.
• This wet/dry handling is a practical engineering screening shortcut, not an unsaturated constitutive model.

12.3 Mesh controls, coupling, and public outputs

The seepage mesh is built with the same shared section mesher used by the deformation workspace. The global target element area can be left on automatic sizing or entered manually, and each polygon may apply a local coarseness factor to refine only that region.

When the engineer enables Use FEM pore pressure in Bishop, the slice base samples pore pressure from the seepage field rather than from the drawn hydrostatic phreatic line. Outside the seepage mesh, or if no seepage result exists, the solver falls back safely to the hydrostatic route.

• Polygon coarseness is a mesh-density control only; it is not a hydraulic soil property.
• The seepage canvas now supports the same high-contrast contour style as deformation: filled contours, isolines, and an in-canvas legend for h, u, |∇h|, |q|, q_x, and q_y.
• The current results panel can plot head, pore pressure, hydraulic gradient, specific discharge magnitude, discharge components, and normal cross-flow along the shared measurement line.
• The line-probe graph is sampled from the solved field and can be copied to the clipboard as distance/value data.
• The seepage workspace is strongest as a 2D screening tool and as a pore-pressure source for the circular Bishop/Spencer solver; it is not yet a full groundwater package.

13.1 Present constitutive route and engineering meaning

The current deformation tool offers two Mohr-Coulomb routes on the shared section mesh: the conservative Stage 1 MC-active reduced-stiffness screen and the Stage 2 exact Mohr-Coulomb elastoplastic return-mapping solve.

Stage 2 is now the default deformation solve for first-pass runs, while Stage 1 remains available as the more conservative hotspot screen. Stage 2 adds stored plastic strain, non-associated flow through ψ, exact Mohr-Coulomb face-edge-apex return logic, exact tension cut-off branches, and an optional plastic geostatic equilibration phase on the shared plane-strain T3 mesh.

The Stage 1 branch is conservative across the monotonic load ramp: once an element exceeds MC during that screening run, it remains on the reduced-stiffness branch for the rest of the increment sequence.

Tension cut-off is diagnostic-only in Stage 1, but constitutive in Stage 2.

For the current Stage 1 screen, the imported CPT layer model also carries a soil-dependent r_shear factor into the deformation material set. Granular soils keep a milder reduction and cohesive soils a stronger one, and the factor remains editable in the deformation material table.

• The default displacement solve uses Stage 2 exact Mohr-Coulomb elastoplasticity, while the conservative Stage 1 MC-active reduced-stiffness screen remains available as a hotspot check.
• Stage 2 stores plastic strain and performs an exact Mohr-Coulomb active-set return with face, edge, apex, and tension-cutoff handling.
• The tool should be read as long-term drained screening, not as an immediate undrained clay check.
• The current implementation uses the Stage 4/5 E_mc value directly as the elastic modulus E.
• The stored dilation angle ψ is active in Stage 2 plastic flow, and plastic geostatic equilibration may be used to equilibrate the initial self-weight state before service loading.

13.2 Load model, supports, and initial stresses

The mechanical load is currently taken from the shared surface load interval on the terrain. In pressure mode the entered q is applied directly; in total-load mode the app converts the total slab load to an equivalent 2D pressure through the out-of-plane length.

The initial effective stress field is built with a geostatic gravity step on the same shared mesh before the external load increment is applied, but the in-situ horizontal confinement is then reset from K_0,nc so elastic ν controls stiffness rather than at-rest stress. Optional seepage coupling only changes the initial pore-pressure field used in that reconstruction.

If plastic geostatic equilibration is selected, the predictor is corrected by a full self-weight equilibrium solve formulated as a predictor-relative constitutive correction rather than a second gravity replay.

• Bottom boundary: u_y = 0. Side boundaries: u_x = 0.
• A practical screening rule is to keep roughly 5B of width on each side of the load and 5B of depth below it; smaller domains will look too stiff.
• Retaining walls are visible in the shared geometry but are not active mechanical elements in the deformation solve.
• If seepage pore pressures are enabled and a seepage result exists, they are sampled into the initial stress field; otherwise the tool falls back to the drawn hydrostatic phreatic line.
• The gravity step lets slopes and wall-supported sections develop initial shear stress before the load increment; if that initialization fails numerically, the solver falls back to the older flat-ground K₀ route and reports a warning.

13.3 Numerical safeguards and public outputs

The deformation solver uses 3-node constant-strain triangles. That choice is robust and fast on the shared mesh, but it comes with the usual stiffness issues for coarse meshes and near-incompressible Poisson ratios.

The public results panel can display settlement, u_x, u_y, |u|, Δσ_yy, initial and final effective / total σ_yy, initial and final effective / total σ_xx, final τ_xy, η_MC, and accumulated plastic strain. In exact tension-cutoff-active zones, η_MC is suppressed because admissibility is governed by the tension surface rather than by the Mohr-Coulomb shear ratio. The deformation overlay can show filled contours, contour lines, and a toggleable legend, and the same shared measurement line can probe those quantities without covering the canvas.

• Poisson's ratio is capped at 0.49 for numerical stability in the plane-strain T3 solver.
• Load-edge refinement is added automatically beneath both ends of the loaded interval because coarse CST/T3 meshes become too stiff there.
• The constitutive update uses the full plane-strain stress₆ state and full principal-stress decomposition.
• Missing or non-positive E_mc is replaced with a pragmatic fallback value of 1000 kPa, with a warning.
• Polygon coarseness remains a local mesh-refinement factor only.
• The line-probe graph can be copied to the clipboard as distance/value data for external checking.

13.4 Arc-length continuation for c-φ safety

Prescribed-strength continuation parameterises the c-φ safety solve by the imposed Σ_Msf and Newton-iterates for displacement. The prescribed control variable becomes ill-conditioned near a limit point: repeated load-step cutbacks below minLoadStep · 4, line-search stalls, an oscillating plastic active set at a fixed Σ_Msf, and a flattening safety curve before a coherent mechanism develops. The WASM CPU path can then switch to a Crisfield–Riks spherical arc-length continuation in which displacement and continuation parameter are unknowns of the same Newton solve.

The corrector is the standard two-solve form (Crisfield 1981): the same tangent stiffness K is solved against −R and against −∂R/∂λ; the linearised constraint ‖W·Δu‖² + α²·Δλ² − Δs² = 0 gives a scalar continuation correction and the displacement correction follows as du = δu_R + dλ · δu_λ. For c-φ safety ∂R/∂λ is computed by finite difference on Σ_Msf against the same committed material state, central where both probes lie in the admissible domain Σ_Msf ≥ 1 and forward-only at Σ_Msf = 1 to respect the strength-reduction floor. Probe assemblies skip tangent construction and never commit material state.

• requestedContinuationMode exposes strength-control (prescribed Σ default), arc-length (forced), and auto (start in strength control, fall back to arc-length on limit-point diagnostics).
• The predictor sign is selected by dot-product against the previous accepted path increment; positive Δλ is forced on the first step within a phase. Under arcLengthAllowPostPeakSafetyPath the predictor may descend (Σ_Msf decreasing with growing displacement) from the second accepted step onward, for diagnostic display only.
• Reported factor of safety remains the highest stable lower-bound Σ_Msf regardless of post-peak descent; post-peak curve points never replace the peak value in the reported FoS.
• Rejected steps restore committed displacement, material-point state, warm-start iterates, and active-set diagnostics exactly.
• Implemented on the WASM CPU path only; the GPU v2 pipeline retains prescribed strength control.

13A.1 Positioning, intended use, and runtime contract

HS targets the three limitations of Mohr-Coulomb that dominate serviceability deformation prediction on real soils: the constant Young modulus, the linear pre-yield response missing the hyperbolic primary curve, and the absence of a primary-loading versus unload-reload distinction. It resolves all three through stress-dependent stiffness, a Duncan–Chang hyperbolic primary curve, and a separate (typically stiffer) unload-reload Young modulus E_ur.

The plugin is implemented entirely in the WASM CPU solver in material_hs.hpp. The JavaScript reference solver explicitly rejects HS analyses at run-start because HS does not register a JS plugin in the material-plugin.js registry. The GPU v2 pipeline does not support HS in this phase; HS analyses run on the WASM CPU backend that is the default for the application.

C-φ safety reduction with HS is supported: c′, tan φ′, tan ψ, and σ_T are divided by Σ_Msf at every active HS Gauss point, while the reference stiffnesses E₅₀^ref, E_oed^ref, E_ur^ref and the exponent m are not scaled, matching PLAXIS safety convention. The HS region constants M_cap, H_cap and sin φ_cv are re-calibrated at every safety trial when φ′ changes via Σ_Msf. Arc-length continuation in the safety phase honours HS through the generic finite-difference Σ_Msf residual derivative.

• Use linear elasticity for regression baselines.
• Use Stage 2 Mohr-Coulomb elastoplasticity for limit-state and safety-driven analyses where the Mohr-Coulomb envelope dominates.
• Use Hardening Soil when serviceability displacements matter and stiffness-versus-stress behaviour controls the engineering output.
• WASM-only. The JavaScript reference path rejects HS; the GPU v2 pipeline rejects HS.
• HSsmall (small-strain stiffness extension, Benz 2007), 3D and axisymmetric paths, anisotropic stiffness, and unsaturated effects are all out of scope in this phase.

13A.2 Stress-dependent reference stiffnesses

Three user-supplied reference stiffnesses pin the elastic and plastic moduli at a reference pressure p_ref (typically 100 kPa). The actual moduli at a Gauss point depend on the current effective stress state through the power law of Schanz et al. (1999), with the minor principal stress σ′₃ governing the triaxial-style stiffnesses and the major principal stress σ′₁ governing the oedometric stiffness.

Typical values of the exponent m are m ≈ 0.5 for cohesionless sands consistent with the Hardin and Drnevich (1972) low-amplitude shear-modulus scaling, m ≈ 1.0 for soft normally-consolidated clays, and m ≈ 0.7–0.9 for stiff clays and silts. The unload-reload Poisson ratio ν_ur is a constant (typically 0.2) that defines the elastic bulk and shear moduli from E_ur via the standard isotropic-elastic relations.

The denominator c′ cos φ′ + σ′ sin φ′ inside the power-law bracket has a numerical floor of 0.5 kPa to prevent the moduli from collapsing when σ′₃ approaches zero. The tension cutoff usually engages before this floor matters, but the floor is the hot-path safeguard.

13A.3 Duncan–Chang hyperbolic primary-loading curve

On primary deviatoric loading at constant confinement, HS reproduces the Duncan and Chang (1970) hyperbolic stress-strain relation. The axial strain in a drained triaxial test follows a hyperbola asymptotic at q = q_a, with q_a larger than the Mohr-Coulomb failure deviator q_f by the inverse of the failure ratio R_f (typically 0.9).

The initial Young modulus E_i is calibrated so that the secant at q = q_f/2 equals E₅₀, giving the closed-form E_i = 2 E₅₀ / (2 − R_f). Below the failure deviator the response is the hyperbolic curve; at q = q_f the Mohr-Coulomb cone takes over as a hard limit and further loading mobilises plastic shear strain γ^p on the cone.

13A.4 Shear (cone) yield surface and mobilised friction

The cone yield function expresses the hyperbolic stress-strain relation as a yield condition with the plastic shear strain γ^p as hardening parameter. After principal decomposition with σ′₁ ≥ σ′₂ ≥ σ′₃ and q = σ′₁ − σ′₃, the cone reads (Schanz et al. 1999, eq. 5):

σ′₂ does not appear in f^s: the cone is a Mohr-Coulomb-aligned surface that depends only on the major and minor principal stresses. The hardening variable γ^p integrates directly from the plastic multiplier of the cone via dγ^p/dλ_s = 1, which keeps the algebra clean.

The mobilised friction angle is the Mohr-Coulomb mobilisation function evaluated at the current stress state. The critical-state friction angle φ_cv is computed once per region at material setup from the inverse Rowe (1962) relation using the user inputs φ′ and ψ. Rowe stress-dilatancy then gives the mobilised dilation angle ψ_m for the cone flow, with two cutoffs: ψ_m = 0 below the critical-state friction angle (no dilation in the contractive regime), and ψ_m = 0 when the void ratio reaches its maximum value (no further dilation beyond e_max).

Plastic flow on the cone is non-associated and uses the Mohr-Coulomb-type plastic potential of Vermeer and de Borst (1984) with the mobilised dilatancy angle. This produces an asymmetric algorithmic tangent that is solved by the existing GMRES dispatch in the deformation solver — no new linear-algebra machinery is needed.

13A.5 Volumetric cap yield surface and K_0,nc calibration

The cap is an ellipse in the meridional plane (p′, q̃) parameterised by the preconsolidation pressure p_p. The tensile shift p_t = c′ cot φ′ translates the origin of the ellipse to the apex of the Mohr-Coulomb cone in compression-positive p′ space. The generalised deviatoric measure q̃ aligns the cap intersection with the Mohr-Coulomb cone in three-dimensional principal-stress space; for triaxial compression with σ′₂ = σ′₃, q̃ reduces to the standard deviator q.

Flow on the cap is associated, so the plastic potential equals the yield function and the cap-mode plastic strain increment is normal to the cap surface. Cap hardening is volumetric and monotonically non-decreasing: the preconsolidation pressure p_p only grows during cap loading, and unloading inside the cap is purely elastic at the unload-reload stiffness.

The cap shape parameter M and the cap-hardening modulus H_cap are pinned by a two-condition calibration along a K_0,nc normal-consolidation loading path: the lateral stress ratio σ₃′/σ₁′ must equal the user input K_0,nc throughout, and the oedometric tangent stiffness at σ₁′ = p_ref must equal E_oed^ref. The implementation solves this two-condition system iteratively at region setup (5–10 single-element evaluations per region), so the hot path never re-computes M or H_cap. The Jaky (1944) default K_0,nc = 1 − sin φ′ is used when the user supplies zero for K_0,nc.

The Brinkgreve (2007) Appendix-A closed form for H_cap as a function of p_p is documented in the feature specification; the current implementation evaluates H_cap once per region at setup and holds it constant during the analysis. The p_p-dependent refinement is a tracked follow-up improvement.

13A.6 Tension cutoff and yield-surface activity regimes

Tension cutoff in HS uses the same surface and the same project-and-check pattern as the exact Mohr-Coulomb route: the tension surface dominates the HS surfaces, and a violation of the tension condition at the elastic trial is projected onto the tension surface before the HS cone or cap correctors run. Tension-controlled states are reported through the active-surface flag with the value 4.

Per Gauss point at the end of a load step, the active surface set is one of five categories: elastic (both HS surfaces inactive and tension admissible), cone-only (cone corrected to f^s = 0, cap admissible), cap-only (cap corrected to f^c = 0, cone admissible), corner (both HS surfaces corrected simultaneously), or tension (tension surface dominant). The corner case is solved by a 2 × 2 Newton in (Δλ_s, Δλ_c) with a closed-form Jacobian; degenerate corner cases fall back to the single-surface returns of the cone-only or cap-only branches.

The implementation hierarchies and failure codes are documented in the feature specification: cone non-convergence returns failure code 101, cap non-convergence 102, corner non-convergence 103, tension-plus-cone-plus-cap apex non-convergence 104, and plane-strain σ_zz inner Newton non-convergence 105. Any of these triggers a load-step cutback through the existing outer-Newton machinery in run_nonlinear_phase.

• Elastic. f^s ≤ 0, f^c ≤ 0, tension admissible. Elastic trial is committed.
• Cone-only. f^s = 0 (corrected). 1-equation Newton in Δλ_s, non-associated flow.
• Cap-only. f^c = 0 (corrected). 1-equation Newton in Δλ_c, associated flow.
• Corner (cone + cap). 2 × 2 Newton in (Δλ_s, Δλ_c) with closed-form Jacobian.
• Tension. Tension surface dominates; HS surfaces re-evaluated at the projected state and corrected if still positive.

13A.7 Return mapping, plane-strain handling, and globalisation

The shipped return mapping is in 6-Voigt with the out-of-plane normal stress σ_zz as an internal degree of freedom that enforces ε_zz = 0 through 1–2 inner Newton iterations per Gauss point per outer Newton step. The plastic potential gradient and the yield-function gradient are evaluated in principal effective stress space and rotated back to the working frame for the stress-update increment.

Convergence tolerances follow the feature specification: the cone return uses |f^s| < 10⁻⁸ · |q_a| and |Δ Δλ_s| < 10⁻¹² with a maximum of 30 inner iterations; the cap return uses |f^c| < 10⁻⁸ · (p_p + p_t)² with the same iteration cap.

The current implementation uses the elastic (E_ur-based) tangent for the global Newton solve rather than the fully consistent algorithmic tangent. The local return mapping certifies per-Gauss-point admissibility, so global convergence in the displacement norm is still well-defined; the elastic-tangent globalisation simply needs more Newton iterations than a consistent tangent would. The consistent tangent (including the σ′₃-dependence of the stress-dependent stiffness and the cap-hardening sensitivity) is tracked as a follow-up improvement.

Predictor projection at seed time. When the K₀ predictor stress state lies outside the HS cap surface (typical for overconsolidated layers with OCR > 1), the predictor is projected onto the cap before the geostatic phase starts. The projection algorithm is the same 1-equation Newton as the cap-only return, identical in structure to the existing predictor projection for the exact Mohr-Coulomb plugin.

13A.8 Material-point state and result envelope

In addition to the common deformation state (effective stress, plastic strain Voigt-6, accumulated equivalent plastic strain), HS stores three new fields per Gauss point: the plastic shear strain γ^p (cone hardening), the preconsolidation pressure p_p (cap hardening), and the total plastic volumetric strain ε_v^p as the signed sum of cone-mode and cap-mode contributions. A 1-byte categorical flag lastActiveSet records the last active surface combination (0 = elastic, 1 = cone, 2 = cap, 3 = corner, 4 = tension).

The wire format extends both input and output to carry the new fields. On input, each region payload carries 12 HS-specific f64 values (E₅₀^ref, E_oed^ref, E_ur^ref, m, ν_ur, p_ref, R_f, K_0,nc, e_init, e_max, OCR, and one alignment slot). On output, each Gauss point carries the existing 284-byte state plus 32 bytes for the HS-specific fields when any region in the model is HS. A model-wide flag hasHsPayload in the result header makes the presence of the HS payload explicit for the JavaScript decoder.

Result overlays. The accumulated equivalent plastic strain and η_MC contour modes work unchanged for HS Gauss points; the η_MC output for HS uses sin φ_m / sin φ′ as the cone mobilisation ratio. New HS-specific contour modes show γ^p, p_p, ε_v^p, and the categorical active-surface flag.

• State. γ^p, p_p, ε_v^p, lastActiveSet stored per Gauss point.
• Active-surface flag. 0 = elastic, 1 = cone, 2 = cap, 3 = corner, 4 = tension.
• Region constants. M_cap, H_cap, sin φ_cv calibrated once per region at material setup.
• Initial state. p_p⁽⁰⁾ = OCR · σ₁′⁽⁰⁾; OCR = 1 places the cap exactly on the in-situ vertical effective stress (normally consolidated), OCR > 1 shifts the cap outward.
• Geostatic seeding. For OCR > 1, the cap is "ahead of" the current state and stays inactive until subsequent loading reaches the preconsolidation level.

13A.9 Parameter cookbook for typical soils

The parameter ranges below follow the standard PLAXIS HS calibration reference of Brinkgreve (2007). They are starting points for preliminary analysis. Site-specific calibration against oedometer, triaxial, or field stiffness data should always supersede generic values. Supporting defaults: ν_ur ≈ 0.2, R_f = 0.9, p_ref = 100 kPa, K_0,nc defaults to 1 − sin φ′ (Jaky 1944) when not entered explicitly, and OCR = 1 for normally consolidated layers.

13A.10 Current implementation limits and follow-up

The HS plugin is production-ready for uniform-load benchmark cases within the convergence envelope characterised in the implementation rollout. The documented limits relative to the published reference formulation are:

• WASM-only deployment. No JavaScript reference path and no GPU support. JS path rejects HS at run-start with a clear error.
• Elastic-tangent globalisation. The global Newton tangent is the E_ur-based elastic tangent rather than the fully consistent algorithmic tangent. Local admissibility per Gauss point is certified; Newton convergence is slower than with a consistent tangent.
• Constant cap-hardening modulus. H_cap is calibrated once per region from the K_0,nc consistency and held constant during the analysis. PLAXIS uses a p_p-dependent H_cap via the Brinkgreve (2007) Appendix-A closed form; this refinement is tracked as a follow-up.
• Narrow convergence envelope at extreme parameters. Very low confinement combined with low cohesion and high friction angle can drive σ′₃ close to zero, where the power-law denominator hits the 0.5 kPa floor. The tension cutoff usually engages before this matters, but very soft surface layers may need finer load-step settings.
• HSsmall (small-strain stiffness extension, Benz 2007), 3D and axisymmetric paths, anisotropic stiffness, and unsaturated effects are all deferred and tracked as separate feature work.