A Bayesian Approach for Inferring Sea Ice Loads

In the 1980s, multiyear sea ice accounted for approximately 50–60% of the total Arctic sea ice area, but in recent years has only accounted for about 15% [1–5]. This dramatic reduction in multiyear ice is part of a significant change in the character of Arctic sea ice, which continues to thin and diminish. As a result, commercial and military activity in the Arctic is likely to increase, necessitating a better understanding of sea ice dynamics and better tools for predicting ice behavior [6–10].

Predictions of ice behavior are typically generated from computational models of the ice dynamics. Both continuum and discrete element methods rely on rheological models to describe the generally nonlinear relationship between stress and strain within the ice. These relationships are often derived for long temporal and spatial scales relevant to climate scale models, where the rheology represents a homogenization of finescale processes such as fracture and ridging. However, on length scales relevant to operations (≈1 km), the ice behavior is more heterogeneous [11–14] and the behavior of individual features, like ridges and leads, becomes more important. In order to better understand the processes that drive ridging and fracture, it is necessary to understand and model the local rheology of the ice and to characterize the ice dynamics over short temporal and spatial scales. To better our understanding, in situ observations of stress and strain in the ice are needed.

Previous field experiments involving in situ stress observations have been spatially localized and have only employed a small number of sensors for a limited amount of time [15–20]. We argue that improving ice rheologies and mechanical process models will require distributed internal stress measurements that continuously monitor the ice stress state. Larger areas of observation and longer monitoring periods are needed to constrain the rheological parameters and to better characterize processes with data that are currently spatially or temporally sparse. Our proposed methodology can help with this by providing a framework that is inexpensive, easy to deploy, and can potentially leverage existing sensing platforms, e.g., buoys that are part of the International Arctic Buoy Programme [21].

Strain is relatively straightforward to measure because the ice deformation can be measured directly. As Ref. [22] points out, however, it is not possible to measure stress directly. Instead, the deformation or strain of another object placed in the ice with known, typically elastic, mechanical behavior must be employed. Using this idea to indirectly measure internal stresses has a long history, starting in earnest with the early works of Refs. [22–24]. Indeed, sensors based on the original design of Ref. [25] are still produced commercially [26]. Our proposed approach uses the same fundamental concept as these early efforts: measuring the deformation of an elastic body frozen in the ice. Although previous systems were special sensors designed specifically for generating point estimates of the strain, we propose to use larger existing buoy platforms and to provide a more comprehensive view of internal ice stress by characterizing the spatially varying pressure field acting between the ice and the buoy. More precisely, we propose a Bayesian inference framework to infer the ice stress applied to the outer wall of a hollow, steel-walled buoy embedded in an ice sheet given a handful of strain observation on the buoy interior. This is made possible by relating the exterior stresses and observed strain with a high fidelity structural finite element model of the buoy. Rigorous approaches to inverse problems, like the one proposed here, are quite powerful and widely used in diverse fields such as medical imaging, geologic exploration, glaciology, and structural health monitoring to estimate quantities that are difficult, expensive, or impossible to measure directly [27–31]. However, we have seen limited use of such techniques in the development of in situ ice sensors.

It is important to note that we do not model the ice itself but instead couple a structural model of the buoy to the ice via spatially distributed traction boundary conditions. Of course, this problem is ill-posed; the traction boundary condition is in theory a function over position, i.e., an infinite dimensional quantity, and therefore cannot be constrained by a finite number of data points. Even after discretizing the traction, we need to infer a parameterized field with many degrees-of-freedom that cannot be completely constrained by the finite number of available strain observations. In fact, there are many traction fields that will match the strain observations equally well. In addition to being ill-posed, observational noise introduces further uncertainty into the problem. Additional knowledge of the traction field (e.g., smoothness or anticipated values) is therefore needed to help overcome the ill-posedness of the problem. Both the ill-posedness and observation noise are naturally accounted for probabilistically in the Bayesian setting employed below.

To assess the utility of our framework, we use observations from a compression test performed at the US Army Cold Regions Research and Engineering Laboratory (CRREL) in the Spring of 2018. Section 2 describes this experiment, the measurement setup, the Finite Element Method (FEM) forward model, and the Bayesian inference framework. In Sec. 3, results are presented for both synthetic data and actual observations from the CRREL experiment. Key patterns in the resulting traction field are related to events witnessed during the experiment. We provide discussion about the strengths and weaknesses with the inference approach in Sec. 4. We conclude with remarks in Sec. 5.

To test the concept of using a buoy instrumented with strain gauges to estimate sea ice internal stresses, we subjected a thin-walled steel buoy to compression within a laboratory grown saline ice sheet. Computationally, we used a CAD representation of the buoy, provided by the University of Washington Applied Physics Laboratory, to develop a structural finite element model of the buoy. This model allowed us to describe the connection between applied traction fields and observed strain and to ultimately define an inverse problem for characterizing the traction field given actual observations. Section 2.2 describes the experimental facility and setup in more detail. Section 2.3 then describes the specifics of our finite element forward model. Finally, Sec. 2.4 formulates a Bayesian inference problem for connecting this model with experimental data in order to characterize the distribution over various loading scenarios.

The ice compression experiment was conducted in the Geophysical Research Facility located at the Cold Regions Research and Engineering Laboratory in Hanover, NH. The facility is composed of a large outdoor tank (dimensions: width 6.7 m, length 18.3 m, depth 2.3 m) and a retractable roof fitted with a refrigeration system designed to facilitate the growth of ice for subsequent mechanical testing. A compressive load was applied to the ice using a hydraulic ram composed of a set of three hydraulic pistons that can together apply a maximum load of 14 MPa.

As shown in Fig. 1, the buoy was constructed from a cylindrical pressure tank with an outer diameter measuring 76.2 cm and a vertical length of 116.8 cm. The pressure vessel was laterally sliced 23.5 cm from the top of the buoy and outfitted with a bolted flange to provide access to the strain measurement system on the interior of the buoy. The buoy was instrumented with 36 strain gauges arranged in pairs that measured vertical and horizontal strain at 18 locations. The strain gauges were placed in three rings around the buoy and vertically spaced 16.5 cm apart, beginning with the top ring of 12 gauges located 34.7 cm from the buoy top (11.2 cm below the bottom of the buoy cap flange), a middle ring located at 51.2 cm from the top (27.7 cm below the bottom of the buoy cap flange), and the bottom ring located at 67.7 cm from the top (44.2 cm below the bottom of the buoy cap flange). The gauges were uniformly distributed radially in 60 deg intervals (see Fig. 1).

To simulate sea ice, the initial tank salinity was set to 27 ppt, which is slightly lower than typical Arctic sea surface salinity (∼29 − 32 ppt) [32]. This salinity level was selected to account for the brine rejection that typically occurs during sea ice formation. The ice growth began on Dec. 19, 2017, and continued until Feb. 12, 2018, reaching a thickness of 23 cm. At this stage of ice growth, the buoy was placed within a square section that was cut into the ice. The ice growing process resumed until reaching the target thickness of 50 cm. During the ice-growing process, the test basin was covered and refrigerated until the experiment took place to promote accelerated ice growth and to limit melting during warm, sunny periods.

The ice compression experiment was conducted on Apr. 11, 2018, which was a partly cloudy day with an outdoor temperature of 9°C. The configuration and load condition resemble a uniaxial compression test illustrated in Fig. 2. The sides of the ice sheet parallel to the loading direction were free of any confinement. The side adjacent to the hydraulic press had a free slip boundary condition, but the opposite side was fixed. For the test, the hydraulic ram load was gradually increased throughout the experiment until reaching a peak value of approximately 14 MPa. The experiment lasted approximately 4 min with strain measurements recorded every 0.5 s.

This linear elastic finite element model was implemented in FEniCS [33] using a tetrahedral mesh that was generated with the cubit software, using a CAD geometry of the buoy used in the experiment (Fig. 1). The CAD files were provided by the University of Washington Applied Physics Laboratory.

We adopt a Bayesian approach to indirectly characterize the ice pressures p_N(x), p_H(x), p_V(x) applied to a defined region on the exterior surface of the buoy from strain gauge measurements on the interior of the buoy. Bayesian inference is a probabilistic framework that allows us to not only estimate the value of the pressures but also rigorously quantify uncertainty stemming from noisy observations and the inability of a few observations to completely constrain the pressures over the entire surface of the buoy. Let Γ_p ⊂ Γ denote the portion of buoy’s exterior surface that lies between 23.5 and 83.5 cm from the top of the buoy. Over this region, we want to infer the pressure function p(x). We start by treating the external pressures applied to the buoy and the observed strains as random variables, denoted by p(x) = [ p_N(x), p_H(x), p_V(x)]^T and ɛ_obs, respectively. The components p_N(x), p_H(x), and p_V(x) denote the stresses in the normal, horizontal, and vertical directions, as defined in Fig. 3.

Our goal is to determine the conditional distribution of the pressures p(x) given the strain observations ɛ_obs. In general, however, p(x) lies in an infinite dimensional function space, which makes inference more challenging. To overcome this, we instead look for pressure functions p(x) that lie within the span of a linear finite element basis; the same basis used above to represent the load vector f. This restriction allows us to characterize the pressure function p(x) with a finite number of coefficients and thus apply standard Bayesian techniques. Ref. [34] provides more information about working directly in the infinite dimensional setting.

We used a set of synthetic strain measurements that resulted from a known load configuration to verify that the inference framework works as expected and provides a reasonable characterization of the applied pressures. The load configuration used for the verification was composed of three rectangular patches that were located with respect to the x-coordinate axis, which is aligned with the θ = 0 and 180 deg direction, as seen in Fig. 4. On the front side, we applied a normal pressure of 4 MPa in a rectangular patch with dimensions 69.9 cm × 20.0 cm and on the back applied normal pressures of 2 MPa in two rectangular patches with dimensions 34.4 cm × 20.0 cm, placed an equal distance on either side of the x-axis. Figure 4 shows the posterior mean of the normal load configurations. The inferred load pattern shows a region of elevated pressure centered at θ = 0 deg with a maximum pressure of 3.47 MPa. The pressures on the rear side of the buoy reach a maximum of 1.65 MPa. The spatial extents of both front and rear regions of high pressure compare well to the known applied load pattern. The spatial pattern of the variances is indicative of the added information that the strain gauges provide and the differences in the length scales for the horizontal and vertical directions in the Gaussian process used to represent the inferred load fields. The variance tends to be high as one moves farther away from a strain gauge with the decay in information being stronger in the vertical direction. Interestingly, the prior lengthscale is actually shorter in the horizontal direction. The slower increase of the posterior covariance in the horizontal direction is therefore an indication that this strain gauge configuration is more informative about horizontal variations in the applied pressures.

Figure 5 provides a more detailed look at the inferred load configuration for horizontal slices at vertical locations of 28.0 and 35.25 cm below the buoy flange. There were aberrant negative loads, but these negative excursions are very small and do not exceed 0.87 MPa. The peak loads for the fore side of the buoy are 3.31 MPa and located at θ = 357 deg and on the aft side of the buoy are 1.56 MPa and located at θ = 184 deg with the load configuration exhibiting a very slight asymmetry. There were small lobes at θ = 100 deg and 260 deg, but these are relatively small in extent and magnitude.

As described in Sec. 2.2, we monitored the buoy’s strain using foil gauges that were affixed to the interior surface of the buoy as an ice sheet was compressively loaded using a hydraulic ram. Figure 6 shows how both the vertically and horizontally aligned strain measured at 18 locations on the buoy varied through the 300 s long experiment. As the ice was initially loaded, the strains quickly rise while the hydraulic ram applied approximately 7 MPa to the ice until the 73 s mark. During this initial period, the vertical strains measured at the locations normal to the direction of loading, θ = 0 and 180 deg, were all tensile as one would expect for a hollow body pinched in the middle. In contrast, the horizontal strains were mostly compressive.

For time period between 73 and 251 s, the hydraulic ram was set to its maximum output of 14 MPa through the duration of the experiment. During this timeframe, we observed a dramatic increase in strain rate for all the locations, and most locations reached a maximum value between 115 and 143 s. There were, however, a few exceptions to this strain trend. For example, the horizontal strains for the θ = 120 and 300 deg locations, which face away and towards the hydraulic ram, respectively, continued to increase after 115 s, but still exhibited a reduction in strain rate when the other locations reached their maximum strain values. In general, the strains were either purely tensile or compressive through the duration of the experiment. However, at the θ = 0 and 240 deg locations, the vertical strains transitioned from compression to tension, or vice versa, with the transitions occurring between 135 and 161 s. In addition, the horizontal gauge at θ = 240 deg transitioned from compression to tension at 124 s.

The bottom horizontal strain gauge at 0 deg stopped collecting data around 10 s, therefore, was removed from the dataset used for the load inversion.

With the Bayesian inference framework introduced in Sec. 2.4, we used the measured strains of the buoy interior wall to infer the loads on the buoy’s exterior surface. We performed this inference on a subregion of the buoy surface, Γ_p, that bounded the region that the ice was expected to be in contact with the buoy. We independently processed each observation to obtain a time series of applied pressures, as shown in Figs. 7 and 8. The inference results are presented as stress fields in the buoy coordinate system, with surface normal direction (positive pointing toward the buoy interior), the horizontal direction (positive pointing clockwise), and the vertical direction (positive pointing up).

For the initial 73 s, the stress levels on the buoy are very low, which is to be expected since the strains were also small during this time interval. When the hydraulic ram output is increased from 7 MPa to 14 MPa, the magnitude of the inferred stresses begin to increase until the horizontal, vertical, and normal stresses reach maxima of 2.7, 1.1, and 20.8 MPa, respectively, near a time of 125 s into the experiment. The posterior mean stress fields at 100 and 125 s exhibit lobes near 0 deg that are consistent with the ice squeezing the buoy at this location and the formation of a crack at 0 deg. The horizontally aligned stress shows a line of zero stress at θ = 0 deg with a positive and negative region on either side of this zero stress line. This stress divergent pattern about the 0 deg line on the buoy is further evidence that a crack had formed at this location. In contrast, the vertical stress pattern shows a zero stress line aligned horizontally with a negative stress (pushing down) above and a positive stress (pushing up) below this line, i.e., tangential stresses that are converging toward this line. The largest vertical stresses are found along the bottom of the buoy’s flange, as shown at time 120 s in Fig. 9. Between 100 and 125 s, the vertical stresses on the flange near 0 deg are positive indicating that ice was pushing up on flange during these times. As Fig. 11 illustrates, this pressure on the flange agrees with visible observations of the buoy being lifted by the flange during the experiment.

Figure 7 provides a more detailed look at the inferred normal stresses along the middle ring of strain gauges (51.2 cm from the top of the buoy) at 100, 120, 160, and 220 s. We see the large magnitude lobes pressure on either side of 0 deg and how the asymmetry increases with time favoring the sextant between 0 deg and 60 deg. There are region of negative stress throughout the experiment that we attribute to an unloading from a prestress in the buoy due to the freezing in process; therefore, we suspect that there is an unknown biases in the strain measurements that we were unable to determine. The maximum inferred pressure in the lobe located in the sextant between 0 deg and 60 deg was 20.5 MPa and the maximum inferred pressure in the lobe located in the sextant located between 320 deg and 0 deg was 13.2 MPa.

As a sanity check, we calculate the posterior predictive mean strain at the observation locations, i.e., the strain ɛ_post, by using the posterior mean pressures as the boundary conditions for the FEM forward model. We compared these strain values against the actual observed strain values, ɛ_obs, as a verification that the posterior loads were matching the observations. The largest difference between the mean posterior predictive and the observed strains was 3.03 × 10⁻¹¹, and occurred at 145 s, which corresponds to shortly after the time of highest observed strain.

The inferred loads were also used to calculate the expected posterior displacement of the buoy surface. The largest displacements were at 175 s, along the 0 deg axis (17 mm), where the buoy deformed inward. In addition to displacement in the direction of loading, the buoy bowed outward along the 90 deg and 270 deg axes (12 mm and 14 mm, respectively) in response to axial compression. The displacement fields align well with the distribution of the normal loads, where highly compressive loads are found in the same areas on the buoy where inward deformation occurred. The displacement results at 200 s are shown in Fig. 10 to illustrate the final deformation state of the buoy. These results help to shape our understanding of the complex loading scenario that caused the observed strains.

To evaluate how well the inference framework worked for this experiment, we qualitatively correlate the inferred load configurations with key events during the physical experiment. Initially, the buoy was frozen into the ice and had full contact between the buoy and the ice through the ice thickness. As the ice was compressively loaded, the buoy acted as an inclusion in the ice and induced localized ice cracking and deformation around the buoy. In some locations, the buoy decoupled from the ice and transferred the load to highly localized regions on the buoy. At approximately 100 s, a large crack formed in front of the buoy along the direction of loading along with two cracks off the backside of the buoy; one along the 180 deg axis and another emanating from the edge of the buoy near 240 deg. These two cracks on the backside separated a small segment of ice from the rest of the ice sheet, as observed in images (c) and (d) of Fig. 11. The normal loads on the front of the buoy presented in Fig. 8 showing compression in two lobes corroborates the presence of a crack in the loading direction. The divergent horizontal loads in the 0 deg direction also indicate that the two sides of this crack were spreading apart as they pressed against the buoy. This is supported by observations that the crack became wider as the experiment progressed.

The ice did not remain planar throughout the experiment and started to buckle at 83 s, with some sections pushing up nearly 1 m above the initial elevation. As it buckled, the ice sheet lifted the buoy unevenly. The uneven lifting action was partly a result of the major crack that formed through the middle (at 0 deg) of the ice sheet. Wlipping between the ice and buoy was also observed until the ice engaged with the buoy flange (see image (b) in Fig 11). The inferred load configurations show large vertical loads on the flange and are highest upward load on the flange at 120 s in Figs. 8 and 9 consistent with the observations. As described previously, the ice load is higher in the sextant between 0 deg and 60 deg matching the visual observations in the experiment where ice exhibiting more engagement with this sextant.

Another key observation from the experiment was that the small separated section of ice on the back edge of the buoy rose with the rest of the ice sheet as it buckled, which is an indicator that the a significant portion of the buoy’s back side, i.e., between 90 deg and 270 deg, was not in contact with ice for a large portion of the test (image (c) in Fig. 11), which may explain why the normal loads were inferred to be either small or even negative (see Fig. 7). Note that by inferring the loads acting on the buoy, we are not explicitly modeling the ice-buoy contact forces. A loss of contact manifests as a decrease in the load from the initial state of the buoy (i.e., negative load in Fig. 7). The posterior probability of a non-negative load could therefore be used to characterize the probability of ice-buoy contact at any location on the buoy. These results and corresponding observations illustrate how our inference framework could qualitatively capture the dynamic and complex loading conditions within the ice sheet that evolved during the experiment.

The good qualitative comparison between the ice behavior and inferred pressures shows promise for our approach. However, there are certainly limitations as well. Certain aspects of the inferred load distribution, such as the high normal loads on the buoy in the direction of loading, are easily explained with observations from the uniaxial compression test. But, other features of the inferred loads are more difficult to justify using the visual observations, such as the negative stress values on the buoy’s backside. Additional experiments with more careful control of the loading conditions will be a better validation exercise for the approach presented in this paper. Another limitation is the small number of observation points (N_obs = 33) relative to the number of degrees-of-freedom (N_dofs = 35,073), which makes the inference problem inherently more challenging because we are trying to predict stress fields with a very sparse dataset. Of course, more observation points would improve results, but may not be feasible. Preliminary tests leveraging Bayesian experimental design techniques have shown promise in optimizing strain gauge locations, but we leave a rigorous study of this approach to future work.

As mentioned previously, we assumed linear elasticity for the forward model to reduce computational complexity, but the buoy did plastically deform breaking this assumption and possibly causing the large inferred pressures. Therefore, a structural model for the buoy that can accommodate plastic deformation should improve the inference results, but it should be noted that a nonlinear model poses added computational challenges since the posterior would no longer be Gaussian and a sampling strategy such as Markov chain Monte Carlo (MCMC) would be required. High-dimensional MCMC for random fields is an active area of research (e.g., Refs. [36,37]) and could be useful for future work in this setting. However, we expect the typical operating regime of deployed buoy systems to be within elastic limits even though this preliminary experiment resulted in plastic deformation.

Being an indirect way of measuring ice stress, our approach will likely be less accurate than existing sensors, such as vibrating wire stress gauges, that are explicitly designed for high resolution stress measurements. However, our approach has the potential to turn any existing ice buoy system into a stress sensor and could therefore lead to much wider temporal and spatial coverage.

We have demonstrated a Bayesian inference framework that infers external load conditions on a buoy using a network of gauges placed on the interior surface of a steel cylindrical buoy. Using the inferred loads as inputs to a structural model of the buoy, we verified the inferred results by comparing the modeled strain against the observed strains. The Bayesian inference approach presented in this paper provides a robust framework for estimating the quantity of interest, in our case the most likely load configuration, but also the associated uncertainty, which helps assess the quality of the inference results.

Distributed observation of internal ice stress within the Arctic ice pack is crucial for improving the understanding of sea ice dynamics. To date, observations of sea ice internal stress are sparse or non-existent and certainly are not part of an ice monitoring system. We believe the framework described here could allow existing buoy systems to be used as a network of stress sensors that can help capture the highly variable stress conditions in sea ice. Continuous, in situ internal stress measurements will allow us to observe how these stresses develop, evolve, and propagate through the ice pack leading to improvements in our understanding of the mechanics behind local processes like ridging and lead formation. This information will be crucial for improving sea ice models and providing operational guidance in the region. The physical character and dynamic behavior of Arctic sea ice will continue to change in response to the warming climate and instrumentation that can provide new insight into the mechanical behavior sea ice will be essential. We believe the Bayesian inference framework presented in this paper is a promising step in that direction.

We would like to thank David Dyer, Jestoni Orcejola, and Dr. Sarah Webster from the University of Washington Applied Physics Laboratory for providing the buoy used in the experiment and the strain gauge data. We would also like to thank Dr. Andrew Davis for taking the time to provide constructive feedback for this manuscript. This study was supported by the Office of Naval Research Arctic and Global Prediction program (code 32) under award N0001418MP00501. This material is based upon work supported by the National Science Foundation under Grant No. ACI-1550487.

There are no conflicts of interest.

Data provided by a third party listed in Acknowledgment.