## Abstract

The Earth’s climate is rapidly changing and some of the most drastic changes can be seen in the Arctic, where sea ice extent has diminished considerably in recent years. As the Arctic climate continues to change, gathering in situ sea ice measurements is increasingly important for understanding the complex evolution of the Arctic ice pack. To date, observations of ice stresses in the Arctic have been spatially and temporally sparse. We propose a measurement framework that would instrument existing sea ice buoys with strain gauges. This measurement framework uses a Bayesian inference approach to infer ice loads acting on the buoy from a set of strain gauge measurements. To test our framework, strain measurements were collected from an experiment where a buoy was frozen into ice that was subsequently compressed to simulate convergent sea ice conditions. A linear elastic finite element model was used to describe the response of the deformable buoy to mechanical loading, allowing us to link the observed strain on the buoy interior to the applied load on the buoy exterior. The approach presented in this paper presents an instrumentation framework that could use existing buoy platforms as in situ sensors of internal stresses in the ice pack.

## 1 Introduction

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% [15]. 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 [610].

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 [1114] 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 [1520]. 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. [2224]. 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 [2731]. 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.

## 2 Methods

### 2.1 Methods Overview.

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.

### 2.2 Physical Experiment.

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).

Fig. 1
Fig. 1
Close modal

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.

Fig. 2
Fig. 2
Close modal

### 2.3 Finite Element Model.

For each strain gauge location in Fig. 1 in the horizontal (i.e., in the x–y plane tangent to the buoy surface) and vertical strain was observed on the interior of the buoy. To relate this to the tractions acting on the exterior of the buoy, we define a model predicting the buoy’s deformation given a particular traction field. This is done by solving the balance of linear momentum equation with the assumption that inertial effects are negligible. This is a reasonable assumption since the loading rate of the ice on the buoy in the experiment is small. With this assumption, the balance of linear momentum is given by
$∇⋅σ(x)+b(x)=0$
(1)
where σ(x) is the Cauchy stress tensor and b(x) is the body force vector on the buoy. Assuming a linear elastic constitutive model, the stress-strain response follows Hooke’s law
$σ(x)=2μεmod(x)+λtr(εmod(x))I$
(2)
where μ and λ are Lamé parameters and ɛmod is the model strain, which is defined as (assuming small strain)
$εmod(x)=12(∇u(x)+(∇u(x))⊤)$
(3)
where u is the displacement. We prescribe a Young’s modulus of 200 GPa and a Poisson’s ratio of 0.3, which correspond to the material properties of the steel buoy used in the experiment. We then derive the Galerkin weak form of Eq. (1) and arrive at the discretized finite element matrix equations
$Kd=f$
(4)
where K is the stiffness matrix, d is the nodal displacement vector, and f is the external nodal force vector acting on the buoy surface Γ. For the purposes of inference, we only need to know predicted values of strain at the physical strain gauge locations on the buoy. Using the finite element basis functions that define the displacement u with the definition of strain in Eq. (3), we are able to write the strain at the observation locations as
$εmod=BK−1f$
(5)
where B is a matrix that contains the derivatives of the finite element shape functions. The size of the B matrix is constructed such that only the shape functions for the finite elements containing the location of the physical strain gauges are included. It is a sparse matrix with 33 rows, one for each strain gauge, and 203, 661 columns, one for each degree-of-freedom in the finite element mesh. The nodal force vector f can be further broken down into
$f=DAp$
(6)
where D is a matrix of shape function products integrated over the external surface Γ, similar to a mass matrix, p is a nodal pressure vector, and A is a matrix that rotates and stamps the pressure vector from the buoy coordinate system, shown in Fig. 3, into the global coordinate system. The matrix A has 203, 661 rows, one for each displacement degree-of-freedom in the mesh, and 35, 073 columns, one for each degree-of-freedom we wish to infer. The pressure vector p is composed of normal, pN, and tangential components in the horizontal and vertical directions, pH and pV, respectively.
Fig. 3
Fig. 3
Close modal
The matrix A transforms the results between Cartesian global coordinates and the buoy coordinates through the expression
$[px(i)py(i)pz(i)]=[An(i)AH(i)AV(i)][pN(i)pH(i)3pV(i)]$
(7)
where
$An=[−n^1(i)−n^2(i)0],AH=[t^H1(i)t^H2(i)0],AV=[00t^V3(i)]$
(8)

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.

### 2.4 Bayesian Inference.

We adopt a Bayesian approach to indirectly characterize the ice pressures pN(x), pH(x), pV(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) = [ pN(x), pH(x), pV(x)]T and ɛobs, respectively. The components pN(x), pH(x), and pV(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.

Let p (without the (x)) denote a vector containing nodal pressures. Our goal is then to characterize the posterior density π(p|ɛobs). Bayes’ rule allows us to write this density as the product of a prior density π(p) and a likelihood function π(ɛobs|p)
$π(p|εobs)∝π(εobs|p)π(p)$
(9)
where the prior density models information known about the pressures before any observations are available, and the likelihood function provides a way of comparing model predictions with the observations. More information is provided on each of these components below.

#### 2.4.1 Prior Distribution.

To obtain a prior distribution over the finite dimensional vector p, we first consider the continuous pressure function p(x), which is spatially varying field. It is natural to probabilistically describe the load function p(x) as a random field. In particular, we describe the prior distribution over the components pN(x), pH(x), and pV(x) with independent Gaussian processes (see e.g., Ref. [35]). A Gaussian process defines a probability distribution over functions and can be interpreted as the infinite-dimensional analog of the more common multivariate Gaussian distribution. Like a multivariate Gaussian distribution, which is completely defined by its mean vector and covariance matrix, a Gaussian process is completely defined by a mean function and a covariance kernel. The covariance kernel describes the correlation between loads at two points x and x′ while the mean function describes the average pressure at a single point. We use μN(x), μH(x), and μV(x) to denote the mean functions for each component and kN(x, x′), kH(x, x′), and kV(x, x′) to denote the corresponding covariance kernels. Together, these functions define Gaussian process distributions over each component
$pN(x)∼GP(μN(x),kN(x,x′))$
(10)
$pH(x)∼GP(μH(x),kH(x,x′))$
(11)
$pV(x)∼GP(μV(x),kV(x,x′))$
(12)
We further assume that all of the covariance kernels take the form
$k(x,x′)=k12(tan−1[x2x1],tan−1[x2′x1′])k3(x3,x3′)$
(13)
where xi denotes component i of the location x, $tan−1[x2x1]$ is the angle around the buoy, k12 is a 1D periodic covariance kernel defining the meridional (horizontal) correlation of p(x), and k3 is a Matern kernel with ν = 3/2 defining the vertical correlation. The standard deviation (i.e., $k(x,x)$) is set to 4.0 MPa for kN and 0.5 MPa for both kH and kV. These standard deviations were chosen to reflect the low probability that pressures would three times these values due to the mechanical properties of the ice. The lengthscale used in all meridional kernels is $π20$ while the lengthscale used in all vertical kernels is 50 cm. These were chosen based on the anticipated smoothness of the pressure fields.
As mentioned above, in order to use these Gaussian process descriptions with the finite element discretization described earlier, we need to discretize the pressure function p(x). To do this with our Gaussian process prior distribution, we evaluate the mean function at every FEM degree-of-freedom in Γp and evaluate the covariance kernel for every pair of locations. The result is a finite dimensional Gaussian distribution that is more convenient to work with than the infinite dimensional Gaussian process. The density over the vector of nodal pressures p is then given by
$π(p)=N(μp,Σp)=1|2πΣp|exp[−12(p−μp)TΣp−1(p−μp)]$
(14)
where
$p=[pNpHpV],μp=[μNμHμV],Σp=[ΣN000ΣH000ΣV]$
(15)
and μN,i = μN(x(i)), ΣN,ij = k(x(i), x(j)), and x(i) denotes the location of mesh node i in Ωp. The Gaussian density π(p) denotes the prior distribution over the external pressures.

#### 2.4.2 Likelihood Function.

The likelihood function describes the distribution of anticipated observations for a particular load p, i.e., what strains are likely to be observed for a particular load? Our finite element model is a relatively accurate representation of the buoy. We would therefore expect the modeled strain ɛmod, given by Eqs. (5) and (6), to be close to the observed strain ɛobs if the true pressures were used in the model. With an accurate model, the difference between ɛmod and ɛobs is then mostly a result of noise in the observations. We model this noise with a Gaussian random variable $z∼N(0,σz2I)$ and assume the additive form
$εobs=εmod+z$
(16)
The likelihood function π(ɛobs|p) is then a Gaussian density over the difference ɛobsɛmod. In particular,
$π(εobs|p)=(2πσz2)−Nz/2exp[−1σz2‖εobs−BK−1DAp‖2]$
(17)
where the finite element model defined in Eqs. (5) and (6) has been introduced to make the dependence on p explicit.

#### 2.4.3 Posterior Distribution.

Multiplying the Gaussian prior in Eq. (14) with the Gaussian likelihood in Eq. (17) yields the Bayesian posterior π(p|ɛobs) ∝ π(ɛobs|p)π(p), which is also Gaussian and defined by a mean vector μp|ɛ and a covariance matrix Σp|ɛ. To simplify notation, consider the prior predictive covariance
$Σε=(BK−1DA)Σp(BK−1DA)T+σz2I$
(18)
and the Kalman gain
$G=Σp(BK−1DA)TΣε−1$
(19)
The posterior mean and covariance are then given by
$μp|ε=μp+G(εobs−(BK−1DA)p)$
(20)
$Σp|ε=Σp−G(BK−1DA)Σp$
(21)
The posterior mean μp|ɛ describes the most likely pressures while the posterior covariance Σp|ɛ characterizes the remaining uncertainty. Note that the block structure in Σp and A enable the posterior mean and covariance to be efficiently constructed without explicitly building or storing the full covariance Σp. The Gaussianity of the prior and likelihood, combined with the use of a linear model, allows us to compute the posterior analytically rather than having to use computationally expensive sampling methods like Markov chain Monte Carlo.

## 3 Results

### 3.1 Inference Verification With Synthetic Strain Measurements.

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.

Fig. 4
Fig. 4
Close modal

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.

Fig. 5
Fig. 5
Close modal

### 3.2 Experimental Strain Measurements.

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.

Fig. 6
Fig. 6
Close modal

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.

### 3.3 Inferred Loads.

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).

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

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.

### 3.4 Posterior Predictive Strains and Displacement.

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−11, 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.

Fig. 10
Fig. 10
Close modal

### 3.5 Visual Observations of Experiment.

Fig. 11
Fig. 11
Close modal

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.

## 4 Discussion

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 (Nobs = 33) relative to the number of degrees-of-freedom (Ndofs = 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.

## 5 Conclusions

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.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

Data provided by a third party listed in Acknowledgment.

## References

1.
Comiso
,
J. C.
,
2012
, “
Large Decadal Decline of the Arctic Multiyear Ice Cover
,”
J. Clim.
,
25
(
4
), pp.
1176
1193
. 10.1175/JCLI-D-11-00113.1
2.
Screen
,
J. A.
, and
Simmonds
,
I.
,
2010
, “
The Central Role of Diminishing Sea Ice in Recent Arctic Temperature Amplification
,”
Nature
,
464
(
7293
), p.
1334
. 10.1038/nature09051
3.
Jeffries
,
M. O.
,
Overland
,
J. E.
, and
Perovich
,
D. K.
,
2013
, “
The Arctic
,”
Phys. Today
,
66
(
10
), p.
35
. 10.1063/PT.3.2147
4.
Carmack
,
E.
,
Polyakov
,
I.
,
,
L.
,
Fer
,
I.
,
Hunke
,
E.
,
Hutchings
,
J.
,
Jackson
,
J.
,
Kelley
,
D.
,
Kwok
,
R.
,
Layton
,
C.
et al.,
2015
, “
Toward Quantifying the Increasing Role of Oceanic Heat in Sea Ice Loss in the New Arctic
,”
Bull. Am. Meteor. Soc.
,
96
(
12
), pp.
2079
2105
. 10.1175/BAMS-D-13-00177.1
5.
Strong
,
C.
, and
Rigor
,
I. G.
,
2013
, “
Arctic Marginal Ice Zone Trending Wider in Summer and Narrower in Winter
,”
Geophys. Res. Lett.
,
40
(
18
), pp.
4864
4868
. 10.1002/grl.50928
6.
Aksenov
,
Y.
,
Popova
,
E. E.
,
Yool
,
A.
,
Nurser
,
A. G.
,
Williams
,
T. D.
,
Bertino
,
L.
, and
Bergh
,
J.
,
2017
, “
On the Future Navigability of Arctic Sea Routes: High-Resolution Projections of the Arctic Ocean and Sea Ice
,”
Marine Policy
,
75
, pp.
300
317
. 10.1016/j.marpol.2015.12.027
7.
Ho
,
J.
,
2010
, “
The Implications of Arctic Sea Ice Decline on Shipping
,”
Marine Policy
,
34
(
3
), pp.
713
715
. 10.1016/j.marpol.2009.10.009
8.
Melia
,
N.
,
Haines
,
K.
, and
Hawkins
,
E.
,
2016
, “
Sea Ice Decline and 21st Century Trans-Arctic Shipping Routes
,”
Geophys. Res. Lett.
,
43
(
18
), pp.
9720
9728
. 10.1002/2016GL069315
9.
Pizzolato
,
L.
,
Howell
,
S. E.
,
Dawson
,
J.
,
Laliberté
,
F.
, and
Copland
,
L.
,
2016
, “
The Influence of Declining Sea Ice on Shipping Activity in the Canadian Arctic
,”
Geophys. Res. Lett.
,
43
(
23
), p.
071489
. 10.1002/2016GL071489
10.
Stroeve
,
J.
,
Markus
,
T.
,
Boisvert
,
L.
,
Miller
,
J.
, and
Barrett
,
A.
,
2014
, “
Changes in Arctic Melt Season and Implications for Sea Ice Loss
,”
Geophys. Res. Lett.
,
41
(
4
), pp.
1216
1225
. 10.1002/2013GL058951
11.
Schreyer
,
H.
,
Sulsky
,
D.
,
Munday
,
L.
,
Coon
,
M.
, and
Kwok
,
R.
,
2006
, “
Elastic-Decohesive Constitutive Model for Sea Ice
,”
J. Geophys. Res.: Oceans
,
111
(
C11
), p.
03334
. 10.1029/2005JC003334
12.
Hutchings
,
J. K.
,
Roberts
,
A.
,
Geiger
,
C. A.
, and
Richter-Menge
,
J.
,
2011
, “
Spatial and Temporal Characterization of Sea-Ice Deformation
,”
Ann. Glaciology
,
52
(
57
), pp.
360
368
. 10.3189/172756411795931769
13.
Marsan
,
D.
,
Stern
,
H.
,
Lindsay
,
R.
, and
Weiss
,
J.
,
2004
, “
Scale Dependence and Localization of the Deformation of Arctic Sea Ice
,”
Phys. Rev. Lett.
,
93
(
17
), p.
178501
. 10.1103/PhysRevLett.93.178501
14.
Rampal
,
P.
,
Weiss
,
J.
,
Marsan
,
D.
,
Lindsay
,
R.
, and
Stern
,
H.
,
2008
, “
Scaling Properties of Sea Ice Deformation From Buoy Dispersion Analysis
,”
J. Geophys. Res.: Oceans
,
113
(
C3
), p.
04143
. 10.1029/2007JC004143
15.
Johnson
,
J.
,
Cox
,
G.
, and
Tucker
,
W.
,
1985
, “
Kadluk Ice Stress Measurement Program
,”
Proceedings of the Eighth International Conference on Port and Ocean Engineering Under Arctic Conditions
, Vol.
1
, pp.
88
100
.
16.
Comfort
,
G.
, and
Ritch
,
R.
,
1990
, Field Measurement of Pack Ice Stresses.
17.
Tucker III
,
W. B.
, and
Perovich
,
D. K.
,
1992
, “
Stress Measurements in Drifting Pack Ice
,”
Cold Regions Sci. Technol.
,
20
(
2
), pp.
119
139
. 10.1016/0165-232X(92)90012-J
18.
Coon
,
M. D.
,
Knoke
,
G. S.
,
Echert
,
D. C.
, and
Stern
,
H. L.
,
1993
, “
Contemporaneous Field Measurements of Pack Ice Stress and Ice Strain Measurements From SAR Imagery
,”
Engineering in Harmony With Ocean. Proceedings
,
, IEEE, pp.
III31
III36
.
19.
Richter-Menge
,
J.
, and
Elder
,
B.
,
1998
, “
Characteristics of Pack Ice Stress in the Alaskan Beaufort Sea
,”
J. Geophys. Res.: Oceans
,
103
(
C10
), pp.
21817
21829
. 10.1029/98JC01261
20.
Hata
,
Y.
, and
Tremblay
,
L.
,
2015
, “
Anisotropic Internal Thermal Stress in Sea Ice From the Canadian Arctic Archipelago
,”
J. Geophys. Res.: Oceans
,
120
(
8
), pp.
5457
5472
. 10.1002/2015JC010819
21.
Zweng
,
M. M.
,
Boyer
,
T. P.
,
Baranova
,
O. K.
,
Reagan
,
J. R.
,
Seidov
,
D.
, and
Smolyar
,
I. V.
,
2018
, “
An Inventory of Arctic Ocean Data in the World Ocean Database
,”
Earth Syst. Sci. Data
,
10
(
1
), p.
677
. 10.5194/essd-10-677-2018
22.
Cox
,
G. F.
, and
Johnson
,
J. B.
,
1983
,
Stress Measurements in Ice
,
USA Cold Regions Research and Engineering Laboratory
,
Hanover, NH
.
23.
Metge
,
M.
,
Strilchuk
,
A.
, and
Trofimenkoff
,
P.
,
1975
, “
On Recording Stresses in Ice
,”
Proceedings of the IAHR Third International Symposium on Ice Problems
,
USA Cold Regions Research and Engineering Laboratory
,
Hanover, New Hampshire
, pp.
459
468
.
24.
Templeton
,
J.
,
1979
, “
On Measurement of Sea Ice Pressures
,”
Proceedings of the Fifth International Conference on Port and Ocean Engineering under Arctic Conditions
,
Trondheim, Norway
, pp.
73
88
.
25.
Johnson
,
J.
, and
Cox
,
G.
,
1980
,
Proceedings of the Workshop on Sea Ice Field Measurements
, Vol.
80
,
Center for Cold Ocean Resources Engineering
,
St. John’s, Newfoundland, Canada
, pp.
193
207
.
26.
2018
,
Instruction Manual: Model 4350BX Biaxial Stressmeter
,
GeoKon
,
Lebanon, New Hampshire
.
27.
Borcea
,
L.
,
2002
, “
Electrical Impedance Tomography
,”
Inverse Probl.
,
18
(
6
), p.
R99
. 10.1088/0266-5611/18/6/201
28.
Schulz
,
R. B.
,
Ale
,
A.
,
Sarantopoulos
,
A.
,
Freyer
,
M.
,
Soehngen
,
E.
,
Zientkowska
,
M.
, and
Ntziachristos
,
V.
,
2010
, “
Hybrid System for Simultaneous Fluorescence and X-Ray Computed Tomography
,”
IEEE Trans. Med. Imaging
,
29
(
2
), pp.
465
473
. 10.1109/TMI.2009.2035310
29.
Morlighem
,
M.
,
Rignot
,
E.
,
Seroussi
,
H.
,
Larour
,
E.
,
Ben Dhia
,
H.
, and
Aubry
,
D.
,
2010
, “
Spatial Patterns of Basal Drag Inferred Using Control Methods From a Full-Stokes and Simpler Models for Pine Island Glacier
,”
West Antarctica, Geophys. Res. Lett.
,
37
(
14
), p.
043853
.
30.
Steiner
,
B.
,
Saenger
,
E. H.
, and
Schmalholz
,
S. M.
,
2008
, “
Time Reverse Modeling of Low-Frequency Microtremors: Application to Hydrocarbon Reservoir Localization
,”
Geophys. Res. Lett.
,
35
(
3
), p.
032097
. 10.1029/2007GL032097
31.
Vanik
,
M. W.
,
Beck
,
J. L.
, and
Au
,
S.
,
2000
, “
Bayesian Probabilistic Approach to Structural Health Monitoring
,”
J. Eng. Mech.
,
126
(
7
), pp.
738
745
. 10.1061/(ASCE)0733-9399(2000)126:7(738)
32.
Brucker
,
L.
,
DIinnat
,
E.
, and
Koenig
,
L.
,
2015
,
Aquarius L3 Weekly Polar-Gridded Brightness Temperature and Sea Surface Salinity
, Version 5,
USA
,
NASA DAAC National Snow and Ice Data Center
.
33.
Alnæs
,
M. S.
,
Blechta
,
J.
,
Hake
,
J.
,
Johansson
,
A.
,
Kehlet
,
B.
,
Logg
,
A.
,
Richardson
,
C.
,
Ring
,
J.
,
Rognes
,
M. E.
, and
Wells
,
G. N.
,
2015
, “
The FEniCS Project Version 1.5
,”
Arch. Numer. Soft.
,
3
(
100
), p.
20553
.
34.
Stuart
,
A. M.
,
2010
, “
Inverse Problems: A Bayesian Perspective
,”
Acta Num.
,
19
, pp.
451
559
. 10.1017/S0962492910000061
35.
Rasmussen
,
C. E.
, and
Williams
,
C. K. I.
,
2005
,
Gaussian Processes for Machine Learning
,
The MIT Press
,
Cambridge, MA
.
36.
Cui
,
T.
,
Law
,
K. J.
, and
Marzouk
,
Y. M.
,
2016
, “
Dimension-Independent Likelihood-Informed MCMC
,”
J. Comput. Phys.
,
304
(
1
), pp.
109
137
. 10.1016/j.jcp.2015.10.008
37.
Cotter
,
S. L.
,
Roberts
,
G. O.
,
Stuart
,
A. M.
, and
White
,
D.
,
2013
, “
MCMC Methods for Functions: Modifying Old Algorithms to Make Them Faster
,”
Stat. Sci.
,
28
, pp.
424
446
. 10.1214/13-STS421