Quantifying Bearing Capacity Loss for Thin-Ice Defects Using a Paired Finite-Element Correction to Gold's Formula

Text copy prepared from the current manuscript source available to ShearLine.

<PARSED TEXT FOR PAGE: 1 / 30>
Quantifying Bearing Capacity Loss for Thin-Ice Defects Using a
Paired Finite-Element Correction to Gold's Formula
Perry Godsea
, Ryley Beddoea
, Shawn Kennyb
, Robert Frederkingc
, Alanna
Carreiraa
a
Civil Engineering, Royal Military College of Canada, 13 General Crerar Cres, Kingston, K7K 7B4,
Ontario, Canada
b
Civil and Environmental Engineering, Carleton University, 1125 Colonel By Drive, Ottawa, K1S 5B6,
Ontario, Canada
 ᶜ National Research Council Canada, Ottawa, K1A 0R6 Canada, Ontario, Canada
Abstract
Ice roads, drill pads, and other engineered floating ice covers in northern Canada depend on reliable
ice-thickness assessment, yet spatially variable ice growth that produces localized thin zones
complicates allowable-load estimation. In practice, localized thin zones in floating freshwater ice
sheets are commonly treated as a minimum or effective ice thickness however, their mechanical
response remains uncertain because it depends on how the thin zone interacts with the flexural
response field of the loaded ice cover. This study examines the effects of circular thin-zone defects
in freshwater ice sheets through 18 laboratory tests comprising uniform and defectve plates with
prepared circular thin zones of varying radius and offset from a localized quasi-static load.
Breakthrough loads, digital-image-correlation-derived vertical deflection profiles, and crack-path
observations are interpreted using reduced-order plate mechanics and paired finite-element post￾processing. Uniform attached-wall tests establish the baseline Gold coefficient, while each defective
plate is compared with an equivalent uniform reference plate to evaluate whether localized thin zones
are better represented as stiffness-transition defects rather than as minimum-thickness measurements
alone. Results show that local thin zones reduced global breakthrough capacity relative to the
uniform-plate baseline. The observed reduction was most strongly ordered by thickness ratio at the
defect, while defect radius and near-edge proximity were more evident in the computed curvature
and relative stress-amplification fields. Peak deflection at breakthrough was not a sufficient indicator
of defect severity such that deflection-profile shape and crack routing provided stronger evidence of
interaction with the active flexural region. A reduced-order Gold-type correction factor is proposed
in which localized defect effects modify the effective coefficient Aeff, rather than replacing nominal
thickness with a minimum value. The proposed framework is intended as a mechanics-informed
interpretation of localized thin-zone effects, not as a calibrated universal design equation.
Keywords: Finite Element; Gold Formula; Localized Thin Zones; Laboratory Tests; Ice, Laboratory
1. Introduction
Ice roads are a critical infrastructure component in northern transportation corridor, providing
seasonal access for fuel and supplies to remote communities and mining operations. Safe operations
require that floating freshwater ice covers sustain transportation loads knowing that ice covers have
spatially variable thickness, changing thermal conditions, and time-dependent material response. In
engineering practice, allowable-load assessments rely primarily on empirical load-thickness
relations of the Gold type, P = Ah2
 to ensure safety (Gold, 1971; Masterson, 2009). However, this
<PARSED TEXT FOR PAGE: 2 / 30>
approach typically uses an over-conservative singular thickness, often taken as the minimum
thickness (hmin) or the effective thickness heff and does not explicitly account for the structural
significance of localized thin zones relative to the load path.
This condition is exposed by modern high-density spatial profiling, such as ground-penetrating radar
(GPR), which identifies localized minima that individual boreholes may fail to capture. While
higher-resolution surveys reveal these defects, their mechanical significance to producing bearing
capacity estimates remains unresolved because yield and ultimate strength is not well agreed upon,
imposing challenges in evaluating structural defects. To investigate the influence of defects on a
plate without establishing a single material strength, this study postulates that the proximity of a thin
region to the active flexural zone that influences the plate's flexural rigidity and hydrostatic support
(Wyman, 1950; Kerr, 1996), is a primary determinant of bearing capacity loss. Consequently,
localized thin zones in this study are treated as discrete spatial stiffness-transition defects embedded
within a floating continuum.
The experimental program performed in this study examines freshwater ice sheets grown in a
controlled cold-room basin at the Royal Military College of Canada, including both nominally
uniform and defective ice sheets that contain localized circular thin zones of varying radii and offset
from a concentrated load. The resulting dataset includes breakthrough loads, deflection profiles
measured by digital image correlation (DIC), and observations of crack-path development. In
engineering practice, professional judgment is often used to infer a single “effective thickness” used
for design when local thickness variability is present. Such an interpretation is useful, but it does not
by itself, describe the influence of localized transition defects on the mechanical response of a
statically loaded floating plate.
This study aims to improve the static, elastic bearing capacity analysis by improving the empirical
nature of Gold’s coefficient (A) through a plate-mechanics finite element interpretation of localized
ice-thickness variability within a Gold-type load-assessment framework. The central premise is that
localized thin zones should not be treated solely as an interpretation of suitable thickness, but as
discrete flexural-rigidity transitions embedded within a floating ice sheet. This study postulates that
structural behaviour of thin zones can be captured through local thickness ratio, defect radius, and
spatial offset relative to the active flexural zone, rather than on a single minimum-thickness value
alone. This framework is consistent with established ice-sheet mechanics, including radial and
circumferential crack development and sensitivity to boundary conditions and loading history
(Sodhi, 1995; Kerr, 1996). 
The objectives of the study are therefore threefold. First, the study quantifies the reduction in
breakthrough load caused by localized circular thin zones relative to uniform-sheet baseline tests.
Second, it assesses the relative importance of local thickness ratio, defect radius, and defect offset
within the tested geometry. Third, it evaluates a reduced-order approach to represent localized thin￾zone effects using Gold’s empirical parameter A. 
2. Background
Floating freshwater ice subject to localized loading conditions are commonly idealized as a thin plate
resting on an elastic foundation. In the classical homogeneous Governing Differential Equation
(GDE), the vertical response of the sheet is governed by the balance between plate bending and
hydrostatic restoring pressure from the displaced water, which may be expressed as:
<PARSED TEXT FOR PAGE: 3 / 30>
where w is vertical deflection, q is applied load intensity, γw is the specific weight of water, and D is the
flexural rigidity of the plate:
where h is the ice thickness, E is the elastic modulus, and is the Poisson ratio. The ideal deflection
of a thin ice sheet under a load is illustrated below in Figure 1.
Figure 1: Thin floating plate loading on an elastic foundation
The analytical plate solutions developed by Wyman (1950), Westergaard (1923), Kerr (1963, 1996),
provide the reference framework for lake ice mechanics. In practice, because the loading is often
highly localized, the governing equation is treated as homogeneous () everywhere except at the point
of application. The localized load is then imposed as a mathematical singularity or a concentrated
area load, allowing the response field to be constrained by measured deflection (Wyman, 1950). This
approach permits the back-inference of an effective flexural rigidity and the estimation of maximum
bending stress from the observed deformation field, even when the physical ice sheet exhibits local
non-homogeneity. 
An immediate implication of this formulation is that stiffness depends on [L3
] dimensional scaling.
However, a localized thin zone in the field is not a geometric minimum but a discrete reduction in
plate thickness, thereby its flexural rigidity. The spatial influence of such a defect is governed by the
characteristic length :
which describes the lateral extent over which deflection, curvature, and bending moments decay
within or proximal to the active flexural zone (USACE, 2002). 
Under downward loading, the sheet develops a flexural deflection bowl with spatially varying
curvature, inducing tensile stresses on the lower surface. Failure is typically associated with the
development of radial and circumferential crack systems (Sodhi, 1995). If a rigidity transition lies
within this critical stress field, it alters the deflection topography and redistributes curvature,
potentially modifying the failure path and breakthrough capacity.
While classical elastic theory provides a primary structural reference, the time-dependent nature of
freshwater ice must be acknowledged. Studies by Mellor (1979), Sinha (1989) and Schulson (1999)
demonstrate that ice exhibits creep, delayed elasticity, and progressive damage under sustained
loading. The observed laboratory breakthrough loads therefore reflect a staged quasi-static loading
history rather than the full time-dependent evolution of creep, damage, or thermal effects. Existing
literature (Sodhi, 1998; Van der Vinne et al., 2013) has largely focused on uniform-thickness ice
covers; however, these studies did not address localized thin-zone defects as the central experimental
variable.
<PARSED TEXT FOR PAGE: 4 / 30>
Accordingly, the gap in literature which motivates the present work is the absence of a mechanics￾based interpretation of localized thin-zone defects. This study addresses this gap by comparing
nominally uniform baselines with sheets containing circular stiffness-transition defects, interpreting
the results in terms of defect proximity and failure-path interaction.
3. Static Ice Plate Mechanics
3.1 Homogeneous Plate Response
The homogeneous floating-plate formulation in Section 2 provides the reference mechanics for
interpreting the laboratory defect tests. In the present section, that framework is used only to define
the quantities required for the paired uniform-versus-defect comparison: the homogeneous reference
response, the local rigidity transition introduced by the defect, and the dimensionless descriptors
used to compare defect severity across tests.
For a concentrated load P acting on an infinite plate, the maximum deflection occurs at the point of
application:
Equation (4) is used as a homogeneous, small-deflection elastic reference solution derived from the
Wyman-type floating-plate formulation in Equation (1). This formulation assumes a linear
hydrostatic restoring response, , and is therefore appropriate only while the deflection remains within
the pre-overwash flotation regime. It does not represent freeboard exhaustion, local submergence,
overwash, or nonlinear changes in water support after the ice surface approaches the water level.
Therefore, neither Equation (1) nor Equation (4) is used to back-calculate material stiffness from
breakthrough deflection, nor are they used as direct predictions of failure deflection. Breakthrough
response may include cracking, local damage, freeboard effects, finite-deflection membrane action,
and geometric nonlinearity that are outside the linear Kirchhoff-Love reference formulation. In this
study, the homogeneous solution is retained to define the reference response, characteristic length,
and flexural scale for comparing uniform and defective plates. Because , thicker plates distribute the
load over a broader flexural domain, whereas local thinning contracts the effective stiffness field and
increases sensitivity to concentrated loading.
3.2 The Mechanics of Gold’s Coefficient
As distributive flexural plate effects are investigated in this study, Gold’s coefficient A is not
interpreted as a purely geometric or purely empirical constant. Rather, it is treated as a lumped
parameter that reflects flexural strength, load concentration, boundary condition, and the plate-scale
response associated with the applied loading geometry. For an area load of radius a, the non￾dimensional load-distribution parameter is:
This parameter measures the radius of the loaded region relative to the characteristic flexural length
scale of the plate. Small values of correspond to strongly concentrated loading, while larger values
correspond to a more distributed load over the flexural domain.
Within the elastic range, the critical surface load that mobilizes a flexural stress in an infinite floating
plate may be expressed in the generalized form first developed by Kerr (1976), and may be written
generally as:
<PARSED TEXT FOR PAGE: 5 / 30>
where Kα is a dimensionless stress-concentration coefficient associated with the axisymmetric plate
solution where its derivation is provided in Appendix C. At flexural failure, σmax = σf therefore,
And by comparing equalities with Gold’s formula, the coefficient A lumps and associates’ failures
loads to stress concentrations as demonstrated below:
Equation (8) demonstrates that the apparent Gold coefficient is not simply a load-thickness
parameter since it accounts for material properties and load-distribution effects, retaining [ML-1T
-2]
dimensions. The same ice thickness may sustain different loads depending on flexural strength, load
radius, boundary condition, and the ratio of the loaded area to the characteristic plate length. This
point is central to the present study because spatial variations in thickness influence both A and σf.
3.3 Curvature-Based Stress Recovery
Concerning the stress-redistributions within the flexural field, the spatially distributed stress state is
recovered from curvature and bending moment. For an axisymmetric plate, the radial and
circumferential bending moments are
and the associated flexural stresses at the plate surfaces are
These expressions are stated here to connect an important distinction in the interpretation of the
experiments. The laboratory measurements, principally breakthrough load and DIC-based deflection
profiles, do not provide stress directly. They provide breakthrough load, a deformation field from
which curvature, moment, where stress may be inferred. For this reason, the interpretation in
Sections 5-7 relies on curvature and stress-based quantities recovered from finite-element analysis,
rather than on deflection magnitude alone. The more relevant connection is how the defect perturbs
curvature, how that perturbation influences recovered flexural stress, and how these effects relate to
breakthrough capacity and the resulting Gold coefficient.
3.4 Localized Thin-Zone Defect as a Rigidity Transition
To address the curvature-driven flexural effects, this study treats the flexural rigidity as a piecewise
rigidity field:
where
Hence, the local stiffness, or the defect in thin plates can be described as:
<PARSED TEXT FOR PAGE: 6 / 30>
thereby defining the defect severity in structural terms. It shows that the defect cannot be
characterized adequately by thickness contrast alone. A local reduction in thickness immediately
produces a cubic reduction in rigidity, and therefore a potentially much larger perturbation in the
flexural response than would be suggested by thickness ratio alone. The defective sheet is therefore
not a homogeneous floating plate with a modified average thickness. It is a variable-rigidity plate
containing a localized stiffness transition. The governing response may be written schematically as
where denotes the Laplacian plate operator modified by spatially varying rigidity. This expression
is used only to indicate the loss of constant-coefficient homogeneity. The full heterogeneous
Kirchhoff-Love operator would require the spatial variation of bending and twisting resultants to be
treated explicitly. Across the defect boundary, displacement and slope continuity must be satisfied
together, with equilibrium of the relevant bending and shear resultants.
A full closed-form solution for Equation (14) is not developed in the present section. That is not
required for the current objective. It is sufficient here to establish that the localized thin zone acts as
a transition defect in rigidity, and that the measured laboratory response must be interpreted as a
perturbation of the homogeneous plate field rather than as an average-thickness problem.
3.5 Defect Geometry & Structural Significance
The effect of the thin-zone defect is governed not only by its stiffness contrast, but also by its
geometry relative to the applied load. Let a denote the radius of the circular load footprint, r the
radius of the thin zone, x the offset distance from the load to the defect, the nominal plate thickness,
and hd the defect thickness. The distance from the load centre to the near edge of the defect is
therefore , while the load edge to defect-edge spacing is . The defect geometry and rigidity contrast
is therefore described by: , and for subsequent stress-amplification analysis, the primary
dimensionless descriptors are:
These quantities distinguish local thickness contrast, flexural-rigidity contrast, defect size relative to
the characteristic flexural length, and defect proximity to the load-centered active flexural region.
<PARSED TEXT FOR PAGE: 7 / 30>
3.6 Basis of Comparison and Scope of Inference
The experimental laboratory program provides both global and local information. Global response
is represented by principally applied loading conditions, and by measuring breakthrough load and
overall deflection response whereas local response is reflected in the DIC profile shape, the position
of the defect relative to the loaded region, and the observed crack-routing behavior at failure. This
distinction is essential, because breakthrough is a global outcome influenced by a local structural
mechanism.
4. Experimental Program
The uniform-sheet tests were used to establish the control group for the laboratory scale testing,
pertaining to the behavior relating to thickness, breakthrough load, and deflection profile for the
attached-wall floating plate. The boundary-condition tests were used to distinguish the effects of
edge restraint, confinement, and sub-ice reaction conditions from the effects of thickness alone, and
thereby to identify the appropriate baseline family for comparison with the defective cases. In one
configuration, a drainage pipe was introduced to relieve confinement effects beneath the sheet within
the basin. The thin-zone tests were then used to examine how a localized rigidity transition modifies
the measured response as a function of defect radius defect offset x, measured from the centre to the
edge, and local thickness contrast within the characteristic length as shown in Figure 2. In this
way, the experiments separated uniform plate behavior from defect-induced perturbations.
4.1 Cold Room Facility and Test Basin
All experiments were conducted in the refrigerated cold room at the Royal Military College of
Canada. The cold room had plan dimensions of 3.6 m by 2.4 m and a height of 2.2 m, with a semi￾hermetic compressor located outside the room and an internal evaporator-fan system used to regulate
the air temperature. As shown in Figure 3, inside the cold room, a wooden test box with internal
dimensions of 2.0 m by 2.0 m and a height of 0.6 m was constructed for the ice-sheet experiments.
The interior walls and bottom of the box were lined with rigid polystyrene foam board to reduce heat
loss and to limit preferential growth along the edges of the sheet. The foam lining was covered with
black plastic to waterproof the box and provide a dark visual background for imaging. The box was
filled with tap water to a depth of approximately 0.20 m. No dopants were added, since the objective
was to study freshwater ice rather than model sea ice.
Figure 2: Various radii of thin zone locations at a distance x cm relative to load area for different size thin zones.
<PARSED TEXT FOR PAGE: 8 / 30>
4.2 Procedure of Ice Specimen Preparation
A repeatable growth procedure was performed to produce freshwater ice sheets of target thickness
subject to control thermal conditions. Prior to freezing, the box was covered with insulating foam
and cloth to reduce surface disturbance from the cold-room fans and to promote initial quiescent
skim formation. To begin a test, the room temperature was reduced to -C C until an initial ice skim
formed, typically after approximately 5 hours. The room temperature was then increased to -C C,
the covers were removed, and the ice was allowed to grow for approximately 20 hours at an average
growth rate of about 0.52 mm/hr. In some tests, a faster growth sequence at - C C for approximately
11 hours was used due to laboratory scheduling constraints. Testing was performed using
temperatures of -C C to limit growth during loading and to test the sheet near to its melting
temperature, consistent with previous freshwater laboratory studies (Gow et al. (1988), Sodhi (1998)
and van Der Vinne et al. (2013)). 
The resulting ice type was not directly measured. No thin-section, universal-stage, microtome, grain￾size, or crystal-fabric measurements were performed as part of the laboratory program. Therefore,
the ice cannot be formally classified as S1 or S2. However, the growth procedure is relevant to
interpreting material variability. Gow (1986) showed that, under quiet (non-turbulent) freshwater
freezing, spray-seeded water produced S2-type columnar ice with predominantly horizontal c-axes,
whereas unseeded or spontaneously nucleated water produced macrocrystalline S1-type ice with
vertical or near-vertical c-axes. Because the present sheets nucleated spontaneously and were not
spray-seeded, the ice is interpreted qualitatively as naturally nucleated, S1-like freshwater
congelation ice. This classification is not used as a measured material input, but it indicates that grain
size, c-axis orientation, and crystal-scale heterogeneity may contribute to variability in breakthrough
load, crack routing, and deflection response. Because no grain-size or crystal-fabric measurements
were collected, possible microstructural effects are treated as unmeasured material variability rather
than as calibrated explanatory variables. This limitation is relevant because naturally nucleated
freshwater ice may develop large grains or preferred crystal orientations that influence flexural
strength, crack initiation, and crack routing. 
Figure 3: Wooden box inside cold room used to grow ice sheets with loading set-up labelled
<PARSED TEXT FOR PAGE: 9 / 30>
4.3 Thin Zone Preparation
Localized thin zones were introduced by placing circular foam inserts on the surface of the ice after
the initial skim had formed. The foam remained in place during the remainder of the growth period,
insulating and delaying local heat transfer beneath it, thereby producing a thinner circular region in
the finished sheet. After the growth period, the foam insert was removed. In some cases, this exposed
open water, which was allowed to refreeze for up to approximately 2 hours before testing. Each thin￾zone test was characterized by the radius of the defect , its offset from the load, and the resulting
local thickness contrast relative to the surrounding sheet. 
The preparation procedure may also have introduced local material differences in addition to
geometric thinning. By insulating the ice surface and delaying local heat extraction, the foam insert
altered the local thermal growth conditions beneath the defect. Where open water was exposed after
insert removal and then refroze, the final defect region may have experienced a distinct local growth
regime under different thermal boundary conditions than the surrounding sheet, potentially
producing differences in grain size or crystallographic orientation. Roos (1966) showed that grain
growth in ice near the melting point can occur over laboratory timescales and that grain-boundary
migration is influenced by both grain shape and crystallographic orientation. Therefore, the prepared
thin zones cannot be assumed to differ from the surrounding sheet only by thickness. Possible local
changes in grain size, crystal orientation, residual flaws, or crack-initiation tendency are treated as
unmeasured material variability. No thin-section, grain-size, or crystal-fabric measurements were
available to separate these effects from the measured thickness field, and no formal ice-type
classification, such as S1 or S2 following Gow (1986), is assigned to either the defect or surrounding
ice. This limitation is carried forward into the interpretation of breakthrough load and crack-path
variability in Section 6.
4.4 Instrumentation
The two principal measured quantities were applied load and ice-sheet deflection. Load was
measured using a Honeywell JP100 load cell with a 0–100 lb capacity and an accuracy 0.15% of
full scale, mounted to the loading arm. As shown in Figure 4 and Figure 5, deflection was measured
primarily using digital image correlation (DIC). 
A Canon DSLR Rebel T6 camera was mounted above the basin and recorded video at 1920 1088-
pixel resolution and 29 frames per second. For post-processing, images were extracted at one frame
per second. Surface targets for DIC tracking consisted of foam blocks marked with random patterns
and attached to the ice using small circular plastic inserts so that target motion reflected ice
deformation rather than slip along the surface. Target spacing was approximately 0.14 m. Linear
Figure 5: Diagram of ice sheet experimental set￾up inside box in cold room
Figure 4:Test T8 Camera field of view (FOV) with instrumentation labelled.
<PARSED TEXT FOR PAGE: 10 / 30>
variable differential transformers (LVDTs) were also used in the original test program as redundant
point measurements and as a check on the DIC-derived deflections, but the distributed DIC
measurements form the primary basis of the present analysis.
The DIC measurements were used primarily to describe the distributed deflection shape and
deformation pattern of each plate, rather than to directly recover pointwise curvature or local peak
stress. This distinction is important because the DIC target spacing was coarse relative to the loading
radius and some defect dimensions, and therefore the DIC field was not interpreted as resolving local
curvature concentrations at the load edge or thin-zone boundary. LVDT measurements provided
fixed-point checks on the timing and magnitude of the DIC-derived vertical deflection. The DIC
measurements were therefore used to interpret global deformation, breakthrough timing, and defect
interaction, while the curvature and flexural-stress fields discussed in Section 6 were recovered from
FEA rather than by differentiating the DIC displacement field. Accordingly, the DIC data support
interpretation of global deformation and defect interaction, while the numerical analysis supports
interpretation of curvature response, stress amplification, and the elastic defect factor.
4.5 Test Procedures
Two boundary-condition issues were considered in the broader laboratory program: sub-ice pressure
condition and edge restraint. To reduce water-pressure buildup beneath the loaded sheet, a PVC
drainage pipe was installed in the corner of the box with its lower portion below water level and
isolated from the main ice sheet by a piece of foam board. 
Ice forming within the pipe was broken prior to testing so that the pipe provided atmospheric pressure
relief during loading. Preliminary open- and closed-pressure tests showed that relieving trapped
water pressure reduced upward edge effects and better approximated the intended attached-wall
floating-plate condition for the main test series. Water-level rise during loading was not measured
independently; the drainage pipe was used only to relieve trapped sub-ice pressure and reduce
confinement effects during the attached-wall tests. A separate free-floating configuration was also
examined by melting the perimeter of the sheet; however, because this represented a different
structural boundary-value problem, it is not included in the principal attached-wall baseline family
used for defect comparison. 
Figure 6: Results from the typical Uniform Ice Sheet where a) Load cell results with corresponding weights applied, and b) Ice sheet
deflection for the applied loads.
Approximately 30 minutes before testing, the DIC targets, camera, and displacement instrumentation
were installed. Load was then applied incrementally through the circular load plate. Initial loading
stages used 0.25 kg increments at approximately 40-second intervals, followed by 1 kg increments
with approximately 60-second waiting periods. In the final stage, 0.5 kg increments were applied
more rapidly until breakthrough occurred. If flooding occurred due to cracking or loss of freeboard,
the final 0.5 kg increments were applied in rapid succession in order to reduce prolonged creep
effects and drive the sheet to breakthrough. The full experimental program comprised uniform-sheet
baseline tests, boundary-condition comparison tests, and localized thin-zone tests. Only tests that
produced clear breakthrough behavior in ice of acceptable quality were retained for the main
analysis.
<PARSED TEXT FOR PAGE: 11 / 30>
5. Experiment Results
The experimental results are presented here in two stages. First, the uniform-sheet tests are used to
define the baseline breakthrough and deflection response of the attached-wall laboratory plate.
Second, the localized thin-zone tests are compared against that baseline to identify the influence of
defect radius, defect offset, and local thickness contrast. 
Observations of crack development and deflection-profile shape are included in this section as
measured or directly observed response features. Their mechanical interpretation is deferred to
Section 6. Plate-Mechanics Post-Processing.
5.1 Response of Uniform Thin Plates
The original laboratory program retained 18 tests after excluding poor-quality tests: six uniform￾thickness tests and twelve thin-zone tests. Of the six uniform tests, one free-floating configuration
was reported separately because significant submergence and over wash made it a distinct boundary￾condition case. The five attached-wall uniform tests retained for baseline comparison comprised two
Figure 7: Maximum deflection profiles of uniform ice sheets. S is sealed pressure and O is for open water pressure
Figure 8: Maximum deflection of DIC targets at breakthrough for thin zone tests with distances of a) 5 cm and b) 16 cm. Vertical
dashed line represents location of thin zone edge from centre.
<PARSED TEXT FOR PAGE: 12 / 30>
sealed-pressure tests and three pressure-relief tests. The 12 thin-zone tests were conducted under
pressure-relief conditions and form the defective dataset used in the paired comparison.
The uniform-sheet results showed a generally positive relation between ice thickness and
breakthrough load. Across the attached-wall tests, thicker sheets sustained larger loads prior to
breakthrough, as expected for a floating plate whose flexural rigidity scale with . The free-floating
test, despite having relatively large thickness, sustained a lower breakthrough load than the attached￾wall tests, indicating the importance of edge restraint in mobilizing load-carrying capacity. In this
respect, the free-floating response should not be regarded as another member of the same baseline
group, but rather as a separate boundary-condition case.
A representative attached-wall uniform test also showed the characteristic progression expected
under staged quasi-static loading. Small initial load increments produced limited measurable
deformation, while later increments generated larger deflections until breakthrough occurred. The
maximum measured deflection occurred beneath the loaded region, and the full DIC profile showed
a bowl-shaped response extending away from the load center. These observations establish the
baseline deformation pattern against which the defective tests are compared.
5.2 Response of Defective Thin Plates
Twelve thin-zone tests were retained for the uniform plate analysis. In these tests, circular defects of
varying radius were placed at prescribed offsets from the load. The resulting breakthrough loads
showed that the presence of a localized thin zone reduced capacity relative to the attached-wall
uniform-sheet baseline, although the magnitude of reduction varied with defect geometry.
When the thin zone was located close to the load, the measured breakthrough loads were generally
lower than those of the more distant defect cases. This trend was most evident when comparing tests
with similar defect size but different offset. By contrast, defect radius did not produce a uniformly
monotonic trend across the full dataset. Rather, the influence of radius appeared to depend on where
the defect was located relative to the loaded region. This indicates that defect proximity is at least as
important as defect size in determining the severity of capacity reduction and suggests that the effect
of a localized thin zone cannot be represented adequately by its thickness or area alone.
The thin-zone tests also showed that maximum recorded deflection at breakthrough did not vary as
strongly as breakthrough load. In other words, different defect configurations could reach similar
peaks and downward displacement while failing at different applied loads. This observation is
important because it indicates that peak deflection alone is not a sufficient indicator of defect
severity. The more informative signal lies in the shape and spatial extent of the measured profile,
together with the relationship between the defect location and the evolving failure path.
5.3 Deflection Profile and Damage Observations
The DIC measurements showed that both uniform and defective sheets developed the expected bowl￾shaped displacement field beneath the load. However, the defective tests differed from the uniform￾sheet baseline in two important respects. First, the radial extent and local curvature shape of the
deflection bowl varied with defect offset. Second, the defective sheets exhibited crack routing
behavior that was not observed in the same manner in the uniform tests.
In the thin-zone tests, radial cracks typically formed first beneath the loaded region and propagated
outward. At breakthrough, one or more circumferential cracks formed and connected the radial crack
system. In many of the defective tests, the circumferential crack passed through or toward the
localized thin zone, indicating that the defect acted as a preferred failure path. This observation is
significant because it suggests that the defect does not simply reduce the average strength of the
sheet but interacts directly with the local stress field and the route by which failure develops.
<PARSED TEXT FOR PAGE: 13 / 30>
Together, the DIC profiles and crack observations indicate that the influence of the thin zone is
expressed more clearly through local deformation geometry and crack routing than through peak
deflection alone. The significance of defect position is therefore not only reflected in breakthrough
load, but also in how the loaded plate redistributes deformation and reaches failure.
5.4 Experimental Interpretation
The preceding observations indicate that the influence of a localized thin zone is not expressed
through breakthrough load alone. The defective sheets did not simply behave as uniform plates of
reduced thickness. Rather, the response depended on how the localized thin region was positioned
relative to the loaded flexural field. This is consistent with the experimental observation that tests
with broadly similar local thickness ratios could exhibit different breakthrough loads and different
crack-routing behavior, while maximum measured deflection at breakthrough did not provide a
unique measure of defect severity.
This distinction is important because the DIC measurements provide a centerline description of a
fundamentally two-dimensional response. The load was applied through a circular contact area, and
the defect was also circular, but the defect was generally offset from the load. The resulting problem
is therefore not axisymmetric except in the uniform reference case. Once the defect is introduced,
the response depends on the spatial relation between the load patch, the thin-zone boundary, and the
surrounding plate. The defect can alter the curvature field locally without producing a proportional
change in peak deflection. In this respect, the measured deflection profile is analogous to a projection
of the plate response rather than a complete description of the stress field.
The crack observations support this interpretation. In several defective tests, radial cracking initiated
near the loaded region and subsequent circumferential cracking developed toward, through, or
adjacent to the thin zone. This behavior suggests that the defect acted as a local modification of the
failure path, rather than as a uniform reduction of global bearing capacity. A purely scalar
interpretation based on , or a weighted mean thickness would not distinguish between a thin zone
located inside the active flexural region and an identical thin zone located outside it. The appropriate
comparison is therefore between the stress field of a defective plate and the stress field of a same￾thickness uniform reference plate.
For this reason, the post-processing analysis is approached as a paired plate analysis problem. For
each defective test, a uniform reference plate of thickness h0 was analyzed under the same load radius
and support condition, and then compared with a corresponding defective plate containing the
measured defect thickness hd, defect radius rd, and defect offset Rd. The purpose of this comparison
is to identify the observational gaps in flexural failure pertaining to the influence of location and size
on how rigidity transitions modify curvature and stress fields.
6. Plate-Mechanics Post-Processing
Using the plate mechanics framework established in section 2 and 3 to investigate the localized 
rigidity transition, the post-processing analysis performs as a paired plate comparison as shown in 
Figure 9. For each defect test, a fictional uniform reference plate of the same nominal thickness, 
load radius, support condition, and applied load was analyzed first. The corresponding defective 
plate is then analyzed using the measured localized thickness field represented as a non-uniform 
flexural-rigidity field. The defect effect was quantified by comparing the recovered stress field in 
the defective plate with the recovered stress field in the equivalent uniform plate.
<PARSED TEXT FOR PAGE: 14 / 30>
This approach isolates stress redistribution caused by the localized rigidity transition while
preserving the same load, spatial extent, and plate geometry. Regions where the defective plate
develops higher recovered stress than the equivalent uniform reference are interpreted as zones of
relative stress amplification. This approach leverages the piecewise function in Equation (11) using
the homogeneous floating plate’s Governing Differential Equation (GDE) described by Equation (1)
to quantify defect behavior as a mesh-dependent local stress amplification and defect factor. The
amplification factor uses the observed stress () and the stress of a thin, uniform infinite plate ()
proposed by Wyman (1950), as a local stress amplification field populated by:
This treats the piecewise stresses at the defect-area transition as the first of three components within
the defect analysis. This quantity identifies locations where the presence of the thin zone increases
or decreases the stress demand relative to the same-thickness uniform plate. Values of Sσ> 1 indicate
local stress amplification, whereas values of Sσ < 1 indicate that the defective plate carries lower
stress locally than the uniform reference plate at that location. The ratio is interpreted only within
the active flexural field of the reference plate. Accordingly, stress-amplification plots are used to
identify the spatial interaction between the defect and the loaded plate, while scalar knockdown
factors are computed from the maximum stress response within the mechanically active region.
 6.1 Selection of an Appropriate Numerical Model
The applied numerical analysis method in principle could be approximated by multiple different
approaches, examining different behaviors. Common applicable methods used to examine flexure
include Finite Element Methods (FEM) (Fitzgerald et al., 2024; Xie, 2022), and Finite Difference
Methods (FDM) (Hindle and Besson, 2023). Due to pseudo-static loading within the prescribed free￾floating flexural extent of the laboratory specimens, both FEM and FDM are at face value,
appropriate approaches to model flexure. However, FDM is ill-suited to the central geometric and
mechanical features of the present problem for two primary reasons.
First, on a structured Cartesian grid, the circular defect boundary must be represented as a stair￾stepped discontinuity. Because flexural rigidity scales as , even modest local errors in the
representation of the transition boundary can produce artificial curvature concentrations.
Furthermore, as demonstrated by Hindle and Besson (2023), standard FDM schemes that rely on
expanding the flexure equation via the product rule are mathematically inconsistent when subjected
Figure 9: Paired finite-element post-processing workflow used to compare an equivalent uniform plate with a localized thin-zone 
defect plate under the same load and support condition.
<PARSED TEXT FOR PAGE: 15 / 30>
to abrupt, non-differentiable changes in rigidity. Such "whole-station" approximations fail to
maintain equilibrium at the transition boundary, potentially leading to “spurious oscillations or
significant errors in the recovered stress field”, and failing to provide correct convergence solutions
for abrupt variations (Cyrus and Fulton, 1966, 1968).
Second, the stress recovery process in FDM relies on numerical differentiation of the computed
deflection field to obtain curvature. While acceptably appropriate for homogeneous plates, this
procedure becomes unreliable near sharp rigidity transitions where grid-imposed discontinuities
contaminate the gradient recovery. For these reasons, FDM was rejected in favor of the FEM.
One advantage of FEM is that it can leverage an unstructured, conformal mesh that represents the
circular defect boundary naturally and allows for local refinement within the active flexural zone.
This reduces geometry-induced artifacts and ensures the piecewise rigidity field is handled with
numerical consistency. Because Kirchhoff-Love plate bending is fourth order in displacement, a
displacement-only conforming formulation requires C¹-continuous interpolation that standard C⁰
triangular elements do not provide. This difficulty is well recognized in thin-plate finite-element
theory, where assumed-displacement Kirchhoff elements require continuity of both transverse
displacement and slope, and mixed formulations are commonly introduced to reduce reliance on
displacement-only stress recovery (Karamanlidis and Atluri, 1984). By applying the mixed finite￾element strategy, the problem was recast using an auxiliary curvature variable, u = −∇ ² w, so that
curvature fields could be obtained within the weak-form solution rather than solely by post￾processing second derivatives of a displacement field as represented graphically in Appendix A.
6.2 Application of Stress Amplification and Elastic Defect Factor
The mixed finite-element post-processing produces a recovered curvature and flexural-stress field
for both the defective plate and the equivalent uniform reference plate. This approach does not
impose curvature continuity across the defect boundary since D(x,y) changes across the thin-zone
edge, the curvature may change sharply while the corresponding bending resultants remain
equilibrated. Consequently, peak recovered curvature and stress near the transition are interpreted as
mesh- and recovery-sensitive local quantities. The stress-amplification field is therefore used as a
paired defect-versus-uniform redistribution index computed with the same mesh family, load
representation, boundary condition, and recovery procedure, rather than as a mesh-independent local
failure stress. 
Applying this post-processing to the FEA field shown in Figure 9 and described by Equation (16)
produces a spatial stress-amplification field that describes the local redistribution of recovered
flexural stress caused by the thin-zone rigidity transition. However, the Gold-type load-thickness
framework introduced in Section 3.2 expresses bearing capacity through the coefficient A in .
Therefore, to connect the spatial FEA result to a Gold-type coefficient modifier, the stress￾amplification field must be condensed to a scalar measure that represents the governing increase in
recovered stress demand. Equation (16) identifies whether the defective plate develops greater or
lower stress demand than the same-thickness uniform reference plate at each point in the domain,
while Equation (17) extracts the maximum stress-amplification factor within the active flexural
region:
The maxima are evaluated within the active flexural region, excluding far-field areas where the
uniform reference stress approaches zero. Normalizing by P expresses the comparison as stress
demand per unit applied load. 
The corresponding elastic defect factor is then defined to translate this maximum stress amplification
into a Gold-type coefficient correction:
<PARSED TEXT FOR PAGE: 16 / 30>
Because is derived from the maximum recovered value of , it inherits the mesh and recovery￾dependence of the stress-amplification field. It is therefore treated as a stable paired elastic defect
factor for the retained mesh family, not as a mesh-independent local material strength or calibrated
design coefficient. Since is applied multiplicatively to the Gold coefficient, any numerical
sensitivity in is carried directly into the defect-modified coefficient , separate from any uncertainty
in the selected baseline coefficient . This mesh sensitivity is evaluated in Section 6.4 before the
factor is used to interpret defect severity.
A value of < 1 indicates that the defective plate develops a higher maximum stress per unit load than
the uniform reference plate, while = 1 indicates that no stress-based reduction is assigned. This scalar
factor is applied in Section 7 to define a defect-modified Gold coefficient, , for evaluating bearing￾capacity loss associated with localized thin-zone defects.
6.3 Numerical Modelling Results
Representative finite-element results are shown for the uniform reference case U1 and the defective
case T7, with the full set of modelled tests provided in Appendix A. The U1 solution provides the
expected axisymmetric homogeneous floating-plate response showing the maximum stress and
curvature auxiliary are concentrated beneath the load and decay radially away from the load centre.
The corresponding amplification field remains essentially uniform because the comparison is made
against the same homogeneous reference state. 
In contrast, the T7 defective solution shows that the stress field remains load-centered in its highest￾magnitude region, but the amplification field is redistributed laterally by the offset thin zone. The
defect does not simply create a uniform stress increase inside the reduced-thickness region, rather,
it perturbs the active flexural field and produces a spatially asymmetric stress-amplification pattern
around the defect boundary. 
Figure 10: a) Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference test
U1showing the expected radially symmetric load-centered response.
Figure 11: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference 
T7, showing latera stress redistribution & transition zone stress amplification.
<PARSED TEXT FOR PAGE: 17 / 30>
This behavior is consistent with the interpretation developed in Sections 6.1- 6.3: the localized thin
zone acts as a rigidity-transition defect whose significance depends on its interaction with the loaded
plate field, not on minimum thickness alone. Because the FEM stress field is used here as a
curvature-based stress proxy rather than a complete fracture criterion, the plots are interpreted
comparatively. Their purpose is to identify and demonstrate the spatial redistribution caused by the
defect and to support the scalar stress-amplification factor used in the reduced-order Gold
interpretation. Full stress, amplification, and curvature-proxy plots for all tests are included in
Appendix A.
6.4 Significance of Stress Amplification and Elastic Defect Factor
The finite-element stress fields provide a spatial interpretation of the defective plates, but a
comparison is also required to assess whether defect severity is controlled primarily by defect size
or by proximity to the loaded region. For this purpose, three normalized descriptors were considered:
the local thickness ratio hd /h0 , the normalized defect radius r/λp , and the normalized near-edge
proximity:
The finite-element factor is the elastic defect factor derived from the computed stress-amplification
field. Its value can be influenced by mesh resolution, interpolation order, and stress-recovery
procedure because it is based on the maximum recovered stress amplification near the thin-zone
boundary. This numerical sensitivity arises because the defect edge is represented as a sharp
transition in , where curvature may change abruptly while the bending resultants remain equilibrated.
In the retained FEM implementation, the transition was represented by assigning and by subdomain
along the conformal defect boundary, rather than by averaging the rigidity across a transition
element. The recovered peak curvature and stress are therefore resolved over the local element scale.
To quantify this effect, was recomputed over the tested mesh-density range. The resulting variation
in was to , with an average variation of approximately . The reported values are therefore treated as
numerically stable defect factors for the applied mesh, but not as mesh-independent local failure
stresses or calibrated design coefficients.
The observed and modelled responses were compared using two separate factors where represented
the measured reduction in breakthrough capacity, while represented the finite-element elastic defect
factor derived from the stress-amplification field. The rank comparison shows that, within the
present test matrix, the observed defect reduction is most strongly ordered by thickness ratio. The
Spearman coefficient between and was , whereas the corresponding values for near-edge proximity
and normalized defect radius were and , respectively. By contrast, the modelled stress factor showed
a stronger geometric sensitivity, with for , for , and for . These results indicate that defect radius
and proximity influence the computed stress redistribution, but they do not independently explain
the observed breakthrough reductions in the available dataset.
This distinction is important for interpretation. The FEM model isolates the elastic plate-mechanics
<PARSED TEXT FOR PAGE: 18 / 30>
effect of the measured geometry, whereas the observed breakthrough load includes ice quality,
growth history, loading sequence, defect preparation, refreezing, crack initiation, and boundary￾condition effects. The rank screening therefore supports a cautious conclusion: thickness contrast is
the clearest observed ordering variable in the present dataset, while defect size and proximity remain
mechanically relevant because they alter the computed stress field. The results do not justify a
calibrated statistical law for offset versus radius, but they support retaining , , and as the appropriate
descriptor set for future defect-sensitive testing.
6.5 Reduced-Order Engineering Gold’s Formula
Because industry-wide safety protocols have integrated the Gold (1971) formula in safety guides, it
is therefore practical to calibrate A0 using the laboratory test results to reinforce inherent risks
associated with A (OIHSA, 2014). The objective is to translate the paired finite element post￾processing results into a form that remains compatible with established load-thickness assessment.
Gold-type relations are retained because they provide a practical empirical foundation to relating
allowable load to ice thickness however, the preceding analysis shows that a localized thin zone
should not be represented only by substituting a minimum thickness into the intact relation. The
defect modifies the stress field of the floating plate, and its influence depends on how the reduced￾rigidity region interacts with the active flexural zone.
In the defective case, the same form is retained, but the coefficient is reduced through the stress￾amplification factor obtained from the plate model. The defect-modified load expression is written
as
where
This preserves the operational simplicity of Gold’s equation while changing the interpretation of the
coefficient. In the intact case, A0 represents the uniform-plate baseline. In the defective case, Aeff
represents a reduced coefficient associated with the computed stress amplification caused by the
localized thin zone. The defect is therefore introduced as a mechanics-informed reduction of intact
static load capacity, not as a practice-based minimum-thickness substitution.
For comparison with the experiments, the observed reduction may be expressed as:
representing the measured reduction of the defective test relative to the same intact load-thickness
baseline. By contrast, is the modelled reduction obtained from the finite-element stress amplification
field. The comparison between and therefore separates the observed capacity reduction from the
elastic stress-redistribution component predicted from the measured defect geometry. Exact
agreement is not expected, because contains experimental effects such as ice quality, growth history,
loading sequence, defect preparation, and crack development, whereas contains only the elastic
plate-mechanics effect of the measured geometry.
To place the laboratory-scale tests in a similar context as the Gold (1971) field data, the load￾thickness plot uses an A-preserving graphical transformation using a scaling parameter s, where:
This transformation preserves the equivalent Gold coefficient because:
<PARSED TEXT FOR PAGE: 19 / 30>
This transformation is purely graphical and -preserving. It is not a physical similitude scaling and is
not interpreted as scaled-model validation.
7. Results
Relative to the Gold-type load-thickness framework introduced by Gold (1971) , the proposed 
stress amplification factor and elastic defect factor are normalized using the Gold A-coefficient￾preserving graphical transformation factor s = 30 to transform the millimetre-scale laboratory data 
into the visual range of field-scale ice-road data from Frederking and Gold (1976). This 
transformation is used only for visualization such that it preserves and therefore preserves the 
equivalent Gold coefficient for each test. The purpose is to compare the laboratory-scale results 
with field-scale data and safe-working-load reference curves without changing the inferred A 
values. The defective test results are shown together with FEM-modified predictions in which the 
stress-amplification-derived elastic defect factor is expressed as an effective coefficient, , as shown
in Figure 12.
The comparison shows that the proposed correction regularizes the defective tests toward the
uniform-plate Gold framework, while still distinguishing observed capacity loss from the elastic
stress-redistribution component predicted by the plate model. The OIHSA (A=3-6) and USACE
(A=10) curves shown in Figure 12 are included only as safe-working-load references and are not
interpreted as failure envelopes. Their purpose is to provide operational context for the Gold-type
load-thickness comparison, while the laboratory points and Gold failure data are used to interpret
breakthrough behaviour.
Table 1: Configuration summary, observed defect factors, FEM-derived elastic defect factors, and Gold coefficients used for the
plotted data in Figure 12.
Test O/S hCentre
(mm)
hdefect
(mm)
X
(cm)
r
(cm)
A
 (kg/cm²)
Pcr
(kN) Bobs Bσ
Figure 12: Gold-type load–thickness context with scaled laboratory results and FEM-modified defect predictions. OIHSA and 
USACE curves are shown as safe-working-load references only and are not interpreted as failure envelopes.
<PARSED TEXT FOR PAGE: 20 / 30>
U1 S 8.50 - - - 11.8 0.0833 - -
U2 S 10.00 - - - 10.7 0.1054 - -
U3 O 10.00 - - - 13.8 0.1354 - -
U4 O 7.50 - - - 17.1 0.0942 - -
U5 O 10.70 - - - 14.5 0.1629 - -
T1 O 9.00 8.95 5 4 13.7 0.0965 0.88 1.00
T2 O 8.40 8.44 5 8 13.6 0.0781 0.82 0.99
T3 O 10.50 10.49 5 16 13.6 0.1262 0.85 0.98
T4 O 10.10 8.98 10 2 13.8 0.1622 1.00 1.00
T5 O 9.50 9.62 10 4 13.7 0.0924 0.76 1.00
T6 O 9.80 9.78 10 8 13.6 0.0924 0.71 0.98
T7 O 16.20 15.90 10 8 13.6 0.2726 0.77 0.99
T8 O 9.60 9.67 10 4 13.8 0.1079 0.86 1.00
T9 O 10.20 10.20 10 16 13.7 0.1508 1.00 0.99
T10 O 11.20 11.15 16 8 13.8 0.1461 0.86 1.00
T11 O 11.70 11.64 16 8 13.8 0.2060 1.00 1.00
T12 O 13.20 13.22 16 15 13.7 0.1854 0.79 0.99
Notes: A is reported to one decimal; βobs and βσ are reported to two decimals. hdefect is the measured local thin-zone thickness used
in the FEM/Gold comparison, while is the nominal reference thickness for the same test. Small differences reflect the measured
defect-thickness field and retained comparison definition.
8. Discussion
The results support the central premise that localized thin zones should not be interpreted only as
minimum-thickness measurements. The uniform baseline tests produced equivalent Gold
coefficients ranging from approximately 10.7 to 17.1 kg/cm², with the attached-wall baseline
centred near the calibrated value used for comparison. The defective tests produced a narrower
equivalent A range of approximately 13.6 to 13.8 kg/cm² after FEM-modified regularization,
indicating that the reduced-order correction brings the defective responses back toward a consistent
Gold-type load–thickness framework. This result is important because it suggests that the Gold form
itself does not need to be discarded; rather, its coefficient must be modified when a localized thin￾zone defect perturbs the plate stress field.
The observed defect factor shows that the localized thin zones reduced measured capacity relative
to the intact baseline. Across the defective tests, ranged from approximately 0.71 to 1.00, meaning
that the observed breakthrough capacity was reduced by as much as about 29% relative to the
baseline Gold prediction. The largest observed reductions occurred in T6, T5, T7, and T12, with
values of approximately 0.71, 0.76, 0.77, and 0.79, respectively. Other tests, such as T4, T9, and
T11, had , indicating that no measurable reduction was assigned under the bounded reduction-factor
definition. This confirms the first objective: localized circular thin zones can reduce breakthrough
capacity, but the magnitude of reduction varies substantially between tests.
The FEM-derived stress factor behaved differently. Across the defective cases, remained close to
unity, ranging approximately from 0.98 to 1.00. This indicates that the elastic plate model predicted
relatively modest stress amplification for the measured defect geometries. The difference between
and is meaningful. It shows that the observed breakthrough reductions contain more than the elastic
stress-redistribution effect alone. Ice quality, defect preparation, refreezing, crack initiation, loading
sequence, and boundary-condition variability likely contributed to the broader spread in . Therefore,
should be interpreted as the elastic plate-mechanics component of defect severity, not as a complete
breakthrough predictor.
The second objective concerned the relative influence of defect offset and defect radius. The rank￾based screening indicates that radius and proximity influence the computed stress field, but they do
not independently explain the observed breakthrough reductions in the available dataset. The
observed reductions were most strongly ordered by thickness ratio , with a Spearman coefficient of
<PARSED TEXT FOR PAGE: 21 / 30>
approximately . By contrast, the observed association with near-edge proximity was negligible, ,
and the association with normalized radius was also weak, . However, the FEM-derived showed
stronger geometric sensitivity, with for near-edge proximity and for normalized defect radius. This
means that radius and offset are mechanically relevant because they affect stress redistribution, but
the current test matrix does not statistically justify that offset alone dominates observed capacity
reduction.
This finding refines the original hypothesis. The results do not support simply that defect offset is
the dominant control on breakthrough load, rather, thickness contrast is demonstrated to be the
clearest observed ordering variable, while defect radius and proximity control how the thin zone
interacts with the computed active flexural field. This is why the appropriate descriptor set remains
as shown in Equation (15), rather than thickness ratio, defect area, radius, or offset alone.
Beyond quantifying the magnitude of capacity reduction, the results also show that peak deflection,
deflection-profile shape, and crack-path observations provide complementary rather than redundant
evidence of defect severity. Peak deflection at breakthrough did not uniquely identify the most severe
defect cases. Several tests reached comparable displacement levels while failing at different loads.
This is consistent with plate mechanics because bending stress is governed by curvature and moment,
not by displacement magnitude alone. The DIC profile shape and crack-routing observations
therefore provide important additional evidence. The tendency for cracks to pass through, toward,
or around the thin-zone transition indicates that the defect affected the failure path and was not
merely a passive local thickness minimum.
A related limitation is that finite-deflection membrane action was not modelled. The Kirchhoff-Love
formulation used in the finite-element post-processing represents a linear bending response and does
not include finite-deflection in-plane stretching or membrane resultants that may develop as
deflection increases relative to ice thickness and freeboard. Therefore, measured breakthrough
deflection is not used to back-calculate elastic stiffness or to validate a nonlinear deflection-to-failure
model. The FEM results are interpreted only as paired elastic stress-redistribution comparisons
between uniform and defective plates under the same modelling assumptions. This distinction is
important near breakthrough, where cracking, freeboard exhaustion, water redistribution, and
membrane action may all contribute to the observed load-deflection response.
The third objective was to evaluate reduced-order representations of localized thin-zone effects. The
proposed FEM-modified Gold expression Equation (20) provides a practical way to retain the Gold
load-thickness form while introducing a mechanics-based defect correction. This is preferable to
replacing h0 with hd , using a weighted average thickness, or treating defect area as the only severity
measure. The FEM correction is modest in the present dataset, but it provides a mechanically
consistent structure where the defect modifies the effective Gold coefficient through the stress field,
not through a simple thickness substitution.
Overall, the results show that localized thin zones are best treated as rigidity-transition defects. The
experimental data demonstrate measurable capacity reductions, the DIC and crack observations
show that defect influence is expressed through deformation shape and failure routing, and the FEM
fields provide a mechanics-based method for converting stress redistribution into an effective Gold
coefficient. The present dataset is not sufficient to produce a universal offset-radius design law, but
it establishes the framework required for one. Future tests should vary , , and independently so that
the elastic stress factor can be calibrated against observed breakthrough reduction with less
experimental coupling.
9. Conclusions
Localized thin zones in floating freshwater ice sheets are better interpreted as rigidity-transition
defects than as simple minimum-thickness measurements. The experiments showed measurable
<PARSED TEXT FOR PAGE: 22 / 30>
reductions in breakthrough capacity, with ranging from approximately 0.71 to 1.00, but the
reductions were not explained by defect radius or offset alone. Within the present dataset, thickness
ratio provided the clearest observed ordering, while radius and near-edge proximity remained
mechanically relevant because they influenced the computed stress redistribution.
Peak deflection at breakthrough was not a sufficient indicator of defect severity. The DIC profile
shape, curvature redistribution, and crack routing provided more useful evidence of how the thin
zone interacted with the active flexural and failure field. The FEM-derived stress factor should
therefore be interpreted as the elastic plate-mechanics component of defect severity, not as a
complete breakthrough predictor.
The proposed reduced-order form, provides a practical bridge between Gold-type load–thickness
assessment and localized plate-defect mechanics. It retains the intact Gold baseline while introducing
the defect through an effective coefficient, , rather than through a minimum-thickness substitution.
Future work should expand the test matrix by varying defect area spatial distribution and geometries,
, , and independently. Additional tests should also place defects at multiple angular positions and
across the plate plan area, rather than along a single line or axis, so that the influence of two￾dimensional defect position and crack-path interaction can be evaluated more completely.
For field application, future research should connect this framework with spatial ice-thickness
mapping methods such as ground-penetrating radar. GPR and related geophysical tools can identify
large-scale thickness anomalies, local depressions, and spatially continuous thin zones along
operational ice roads. The mechanics framework developed here suggests that such anomalies should
not be screened by minimum thickness alone. Instead, mapped thin zones should be assessed
according to their thickness contrast, spatial extent, and proximity to anticipated vehicle load paths.
In this way, high-resolution thickness data could be converted into defect-sensitive load-reduction
maps for operational decision-making.
Acknowledgement
The authors gratefully acknowledge the Royal Military College of Canada for supporting the laboratory
testing program, staff, students and personnel who contributed towards the experimental program and data
collection that formed the basis of this study. The authors acknowledges mentorship, technical discussions,
and research training received at Carleton University, which contributed to the development of the plate￾mechanics interpretation and freshwater ice-engineering context used in this manuscript. The authors also
acknowledge the National Research Council Canada for its longstanding documentation, testing, and
reporting on freshwater ice engineering, which provided valuable background context for this work, and
support through financial support of the Natural Sciences and Engineering Research Council of Canada
through NSERC Alliance Grant ALLRP 575246-22.
 References
[1] Aalto University, “Aalto Ice Tank,” 2018. [Online]. Available: https://www.aalto.fi/research￾andlearning-infrastructures/aalto-ice-tank.
[2] A. J. Gow, “Orientation textures in ice sheets of quietly frozen lakes,” Journal of Crystal Growth,
vol. 74, pp. 247–258, 1986.
[3] A. J. Gow, H. T. Ueda, J. W. Govoni, and J. Kalafut, "Temperature and structure dependence of
the flexural strength and modulus of freshwater model ice," Cold Regions Research and
Engineering Laboratory, U.S. Army Corps of Engineers, Report 88-6, 1988.
[4] M. Aly, R. Taylor, E. B. Dudley, and I. Turnbull, “Scale effect in freshwater ice flexural
strength,” in Proc. OMAE2018-77314, 2018, 9 p.
<PARSED TEXT FOR PAGE: 23 / 30>
[5] Standard Guide for Calibration of Linear Displacement Sensor Systems Used to Measure
Micromotion, ASTM Standard F2537-06, 2011.
[6] Z. P. Bažant, “Scaling of sea ice fracture – Part I: Vertical penetration,” J. Appl. Mech., vol. 69,
pp. 11-18, 2002.
[7] Z. P. Bažant and Y.-N. Li, “Penetration fracture of sea ice plate: Simplified analysis and size
effect,” J. Eng. Mech., vol. 120, no. 6, pp. 1304-1321, 1994.
[8] S. Beltaos, “Bearing capacity of floating ice covers: theory versus fact,” in Proc. 11th Workshop
on River Ice: River Ice Processes Within a Changing Environment, May 2001, pp. 14-16.
[9] BMT Fleet Technology, “The Rideau Canal Skateway in-situ and laboratory ice property tests
for determining the bearing capacity of the ice,” Report 6441.FR, Prepared for the National
Capital Commission, 2010.
[10] D. M. Cole, “The microstructure of ice and its influence on mechanical properties,” Eng. Fract.
Mech., vol. 68, pp. 1797-1822, 2001.
[11] J. P. Dempsey, “Fracture toughness of ice,” in Proc. IUTAM-IAHR Ice-Structure Interaction
Symp., Berlin, 1991, pp. 109–145.
[12] J. P. Dempsey, S. J. DeFranco, R. M. Adamson, and S. V. Mulmule, “Scale effects on the in￾situ tensile strength and fracture of ice. Part I: Large grained freshwater ice at Spray Lakes,
Reservoir, Alberta,” Int. J. Fract., vol. 95, pp. 325–345, 1999.
[13] S. H. Faria, I. Weikusat, and N. Azuma, “The microstructure of polar ice. Part II: State of the
art,” J. Struct. Geol., vol. 61, pp. 21-49, 2014.
[14] C. Fox, “A scaling law for the flexural motion of floating ice,” in Proc. IUTAM Symp. Scaling
Laws in Ice Mechanics and Ice Dynamics, 2001, pp. 135-148.
[15] G. E. Frankenstein, “Load test data for lake ice sheets,” U.S. Army Cold Reg. Res. Eng. Lab,
Hanover, NH, Tech. Rep. 89, 1963.
[16] G. E. Frankenstein, “Strength of ice sheets,” in Proc. Conf. Ice Pressures Against Structures,
Laval Univ., Quebec, Nov. 10-11, 1966, NRCC No. 9851, pp. 79-87.
[17] R. M. W. Frederking and L. W. Gold, “The bearing capacity of ice covers under static loads,”
Can. J. Civ. Eng., vol. 3, pp. 288-293, 1976.
[18] L. W. Gold, “Use of ice covers for transportation,” Can. Geotech. J., vol. 8, no. 170, pp. 170-
181, 1971.
[19] L. W. Gold, “Engineering properties of fresh-water ice,” J. Glaciology, vol. 19, no. 81, pp. 197-
212, 1977.
[20] K. M. Gow, H. T. Ueda, and J. S. Richter, “Flexural strength of freshwater ice,” U.S. Army
Cold Reg. Res. Eng. Lab, Hanover, NH, 1988.
[21] M. Hetényi, Beams on Elastic Foundation: Theory with Applications in the Fields of Civil and
Mechanical Engineering. Ann Arbor, MI, USA: Univ. of Michigan Press, 1946.
[22] D. Karamanlidis and S. N. Atluri, “Mixed finite element models for plate bending analysis:
Theory,” Computers & Structures, vol. 19, no. 3, pp. 431–445, 1984.
[23] A. D. Kerr, “The bearing capacity of floating ice plates subjected to static, moving and
oscillatory loads,” J. Glaciology, vol. 17, no. 76, pp. 229-268, 1976.
<PARSED TEXT FOR PAGE: 24 / 30>
[24] A. D. Kerr, “Bearing capacity of floating ice covers subjected to static, moving, and oscillatory
loads,” Appl. Mech. Rev., vol. 49, no. 11, pp. 463-476, 1996.
[25] M. Leppäranta, Freezing of Lakes and the Evolution of their Ice Cover. Berlin, Germany:
Springer-Verlag, 2015.
[26] D. M. Masterson, “State of the art of ice bearing capacity and ice construction,” Cold Reg. Sci.
Technol., vol. 58, pp. 99-112, 2009.
[27] M. Mellor, “Mechanical properties of polycrystalline ice,” in Physics and Mechanics of Ice,
Berlin: Springer, 1979, pp. 217–245.
[28] OIHSA, “Best Practices for Building and Working Safely on Ice Covers in Ontario,” Ontario
Infrastructure Health and Safety Association, 2014.
[29] E. M. Schulson, “The brittle-to-ductile transition in ice under compression,” J. Phys. Colloq.,
vol. 48, pp. 413-420, 1987.
[30] N. K. Sinha, “Rheology of column-grained ice,” Exp. Mech., vol. 18, no. 12, pp. 464-470, 1978
[31] D. S. Sodhi, “Breakthrough loads of floating ice sheets,” J. Cold Reg. Eng., vol. 9, no. 1, pp. 4–
20, 1995.
[32] D. S. Sodhi, “Deflection and cracking of floating ice sheets,” Cold Reg. Sci. Technol., vol. 27,
no. 3, pp. 175-190, 1998.
[33] USACE, “Engineering and Design: Ice Engineering,” EM 1110-2-1612, U.S. Army Corps of
Engineers, Hanover, NH, 2002.
[34] USACE, “Safety on Floating Ice Sheets,” U.S. Army Corps of Engineers, Hanover, NH, 2003.
[35] R. Van der Vinne, D. S. Sodhi, and J. P. Dempsey, “Experimental study on the breakthrough of
floating ice sheets,” in Proc. 22nd Int. Conf. Port Ocean Eng. Arctic Cond., Espoo, Finland,
2013.
[36] H. M. Westergaard, “Om beregning af plader paa elastisk underlag med særligt henblik paa
kørebaner i betonveje,” Ingenioren, vol. 32, pp. 513-524, 1923.
[37] M. Wyman, “Deflections of an infinite plate,” Can. J. Res., vol. 28a, no. 3, pp. 293–302, 1950.
[38] D. Hindle and O. Besson, "A corrected finite-difference scheme for the flexure equation with
abrupt changes in coefficient," Solid Earth, vol. 14, no. 1, pp. 197–212, Mar. 2023, doi:
10.5194/se-14-197-2023.
[39] Xie, Q. (2022). A study of numerical simulation on extreme flexural stress of ice plate. Material
Science, Engineering and Applications, 2(1), 9–14. https://doi.org/10.21595/msea.2022.22411
<PARSED TEXT FOR PAGE: 25 / 30>
[40] Fitzgerald, A., & Janse van Rensburg, W. (2023). Limitations of Gold’s formula for predicting
ice thickness requirements for heavy equipment. Canadian Geotechnical Journal, 61(1), 183–
188. https://doi.org/10.1139/cgj-2022-0464
Appendix A
Figure A.1 a) Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference 
test U1 showing the expected radially symmetric load-centered response.
Figure A.2: a) Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference 
test U2 showing the expected radially symmetric load-centered response.
Figure A.3: a) Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference 
test U3 showing the expected radially symmetric load-centered response.
Figure A.5: Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference test
U5 showing the expected radially symmetric load-centered response.
Figure A.4: a) Representative finite-element stress b) stress-amplification and c) curvature-proxy fields for the uniform reference 
test U4 showing the expected radially symmetric load-centered response.
Figure A.6: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference, 
showing latera stress redistribution & transition zone stress amplification.
Figure A.7: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference 
T2, showing latera stress redistribution & transition zone stress amplification.
Figure A.8:
T3, showing latera stress redistribution & transition zone stress amplification.
<PARSED TEXT FOR PAGE: 26 / 30>
Figure A.9: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference 
T4, showing latera stress redistribution & transition zone stress amplification.
Figure A.10: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T5, showing latera stress redistribution & transition zone stress amplification.
Figure A.11: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T6, showing lateral stress redistribution & transition zone stress amplification.
Figure A.12: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T7, showing latera stress redistribution & transition zone stress amplification.
<PARSED TEXT FOR PAGE: 27 / 30>
Figure A.13: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T8, showing latera stress redistribution & transition zone stress amplification.
Figure A.14: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T9, showing latera stress redistribution & transition zone stress amplification.
Figure A.15: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T10, showing latera stress redistribution & transition zone stress amplification.
Figure A.16: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T11, showing latera stress redistribution & transition zone stress amplification.
<PARSED TEXT FOR PAGE: 28 / 30>
Appendix B
Figure A.17: a) Representative finite-element stress b) stress-amplification, and c) curvature-proxy fields for the uniform reference
T12, showing latera stress redistribution & transition zone stress amplification.
Figure B.1: Gold-type load–thickness context with scaled laboratory results and FEM-modified defect predictions.
Figure B.2: Load P vs laboratory ice thickness showcasing the observations from the uniform and defect tests as well as a calibrated 
line of best fit, representing A Figure B.3: Uniform plate results compared to the 18-data point fitted line A 0 = Aeff = 13.81 0 = 13.81
<PARSED TEXT FOR PAGE: 29 / 30>
Appendix C
C.1 Freeboard and Incipient Submergence
The homogeneous floating-plate formulation used in Equation (1) assumes a linear hydrostatic
foundation. In this formulation, the upward restoring pressure is proportional to downward
deflection:
where is the hydrostatic restoring pressure, is the unit weight of water, and is the downward
deflection. This approximation is appropriate while the ice surface remains within the pre-overwash
flotation regime. It does not represent freeboard exhaustion, local submergence, overwash, or
nonlinear changes in water support after the ice surface approaches the water level.
For a floating ice cover, the available freeboard is defined as:
where is ice thickness and and are ice and water density, respectively. Incipient submergence
occurs when the maximum central deflection equals the available freeboard, i.e., when w0 = f. For
the homogeneous plate reference,
Where
Incipient freeboard exhaustion is evaluated by setting the homogeneous reference deflection in
Equation (C.3) equal to the available freeboard in Equation (C.2), and substituting the definition of
from Equation (C.4). Solving for gives:
For GPa, , kg/m3
, kg/ m3 and N/m3
, Equation (C.5) reduces numerically to:
with hi in metres and P in Newtons. 
Equation (C.6) is retained only as a flotation-regime indicator. It defines the hydrostatic limit at
which freeboard is exhausted. It does not define flexural breakthrough, which is governed by stress,
curvature, and fracture development rather than flotation alone.
C.2 Flexural Stress and Load-Distribution Scale
The more important reference for the present paper is the flexural stress state associated with
localized loading. The stress transfer mechanism of the floating plate depends not only on load
magnitude, but on how concentrated the load is relative to the plate length scale. For an area load of
radius, a, the non-dimensional load-distribution parameter is defined as
<PARSED TEXT FOR PAGE: 30 / 30>
This parameter measures the radius of the loaded region relative to the characteristic flexural length
scale of the plate. Small values of correspond to strongly concentrated loading; larger values
correspond to more distributed loading.
Within the elastic range, the critical surface load that mobilizes a flexural stress in an infinite floating
plate may be expressed in the generalized form first developed by Kerr (1976):
where is the derivative of the Kelvin-function term associated with the axisymmetric infinite-plate
solution. Rearrangement gives the corresponding stress-concentration form:
Combining these expressions, the critical load may then be written as:
In comparison with the empirical Gold form:
It may then describe the stress concentration as an equality that is integrated in Gold’s empirical
parameter A, and expressed mechanical interpretation:
The equations shows that the apparent Gold coefficient (A) is not, in plate-mechanics terms, a purely
geometric or purely empirical quantity. Rather, it is a lumped parameter reflecting both the flexural
strength of the ice and the non-dimensional load-distribution effect embedded in .
∗
Corresponding author. – Perry Godse, M.Sc., EIT
August 28, 2026