A Bayesian Approach for Inferring Sea Ice Loads
Matthew Parno, Taylor Hodgdon, Brendan West, Devin O'Connor, Arnold Song
MMatthew Parno US Army Engineer Research andDevelopment Center,Hanover, NH 03755 USAemail: [email protected]
Taylor Hodgdon
US Army Engineer Research andDevelopment Center,Hanover, NH 03755 [email protected]
Brendan West
US Army Engineer Research andDevelopment Center,Hanover, NH 03755 USAemail: [email protected]
Devin O’Connor
US Army Engineer Research andDevelopment Center,Hanover, NH 03755 USAemail: [email protected]
Arnold Song
Dartmouth College Hanover, NH 03755 USAemail: [email protected]
A Bayesian Approach for InferringSea Ice Loads
The Earth’s climate is rapidly changing and some of the most drastic changes can beseen 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 in-creasingly important for understanding the complex evolution of the Arctic ice pack. Todate, 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 withstrain gauges. This measurement framework uses a Bayesian inference approach to inferice loads acting on the buoy from a set of strain gauge measurements. To test our frame-work, strain measurements were collected from an experiment where a buoy was frozeninto ice that was subsequently compressed to simulate convergent sea ice conditions. Alinear elastic finite element model was used to describe the response of the deformablebuoy to mechanical loading, allowing us to link the observed strain on the buoy interiorto the applied load on the buoy exterior. The approach presented in this paper presentsan instrumentation framework that could use existing buoy platforms as in situ sensors ofinternal stresses in the ice pack.Keywords: Sea Ice, Stress Measurement, Bayesian Inference, Model-based Sensing
In the 1980s, multiyear sea ice accounted for approximately50-60% of the total Arctic sea ice area, but in recent years hasonly accounted for about 15% [1–5]. This dramatic reduction inmultiyear ice is part of a significant change in the character ofArctic 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 bettertools for predicting ice behavior [6–10].Predictions of ice behavior are typically generated from compu-tational models of the ice dynamics. Both continuum and discreteelement methods rely on rheological models to describe the gen-erally nonlinear relationship between stress and strain within theice. These relationships are often derived for long temporal andspatial scales relevant to climate scale models, where the rheologyrepresents a homogenization of finescale processes like fractureand ridging. However, on length scales relevant to operations ( ≈ Corresponding AuthorFebruary 18, 2021 rently spatially or temporally sparse. Our proposed methodologycan help with this by providing a framework that is inexpensive,easy to deploy, and can potentially leverage existing sensing plat-forms, e.g., buoys that are part of the International Arctic BuoyProgramme [21].Strain is relatively straightforward to measure because the icedeformation can be measured directly. As [22] points out however,it is not possible to measure stress directly. Instead, the defor-mation or strain of another object placed in the ice with known,typically elastic, mechanical behavior must be employed. Usingthis idea to indirectly measure internal stresses has a long his-tory, starting in earnest with the early works of [22–24]. Indeed,sensors based on the original design of [25] are still producedcommercially [26]. Our proposed approach uses the same funda-mental concept as these early efforts: measuring the deformationof an elastic body frozen in the ice. However, whereas as previoussystems were special sensors designed specifically for generatingpoint estimates of the strain, we propose to use larger existing buoyplatforms and to provide a more comprehensive view of internalice stress by characterizing the spatially varying pressure field act-ing between the ice and the buoy. More precisely, we propose aBayesian inference framework to infer the ice stress applied to theouter wall of a hollow, steel-walled buoy embedded in an ice sheetgiven a handful of strain observation on the buoy interior. This ismade possible by relating the exterior stresses and observed strainwith a high fidelity structural finite element model of the buoy.Rigorous approaches to inverse problems, like the one proposedhere, are quite powerful and widely used in diverse fields such asmedical imaging, geologic exploration, glaciology, and structuralhealth monitoring to estimate quantities that are difficult, expen-sive, or impossible to measure directly [27–31]. However, we haveseen limited use of such techniques in the development of in situice sensors.It is important to note that we do not model the ice itself butinstead couple a structural model of the buoy to the ice via spa-
Journal of Applied Mechanics
PREPRINT FOR REVIEW / 1 a r X i v : . [ c s . C E ] F e b ially distributed traction boundary conditions. Of course, thisproblem is ill-posed; the traction boundary condition is in the-ory a function over position, i.e., an infinite dimensional quantity,and therefore cannot be constrained by a finite number of datapoints. Even after discretizing the traction, we need to infer aparameterized field with many degrees of freedom that cannot becompletely constrained by the finite number of available strain ob-servations. In fact, there are many traction fields that will matchthe strain observations equally well. In addition to being ill-posed,observational noise introduces further uncertainty into the prob-lem. Additional knowledge of the traction field (e.g., smoothnessor anticipated values) is therefore needed to help overcome the ill-posedness of the problem. Both the ill-posedness and observationnoise are naturally accounted for probabilistically in the Bayesiansetting employed below.To assess the utility of our framework, we use observationsfrom a compression test performed at the US Army Cold RegionsResearch and Engineering Laboratory (CRREL) in the Springof 2018. Section 2 describes this experiment, the measurementsetup, the Finite Element Method (FEM) forward model, and theBayesian inference framework. In Section 3 results are presentedfor both synthetic data and actual observations from the CRRELexperiment. Key patterns in the resulting traction field are relatedto events witnessed during the experiment. We provide discussionabout the strengths and weaknesses with the inference approach inSection 4. We conclude with remarks in Section 5. To test the concept of using a buoyinstrumented with strain gauges to estimate sea ice internalstresses, we subjected a thin-walled steel buoy to compressionwithin a laboratory grown saline ice sheet. Computationally, weused a CAD representation of the buoy, provided by the Universityof Washington Applied Physics Laboratory, to develop a struc-tural finite element model of the buoy. This model allowed usto describe the connection between applied traction fields and ob-served strain and to ultimately define an inverse problem for char-acterizing the traction field given actual observations. Section 2.2describes the experimental facility and setup in more detail. Sec-tion 2.3 then describes the specifics of our finite element forwardmodel. Finally, Section 2.4 formulates a Bayesian inference prob-lem for connecting this model with experimental data in order tocharacterize the distribution over various loading scenarios.
The ice compression experimentwas conducted in the Geophysical Research Facility located at theCold 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 rooffitted with a refrigeration system designed to facilitate the growthof ice for subsequent mechanical testing. A compressive load wasapplied to the ice using a hydraulic ram composed of a set of threehydraulic pistons that can together apply a maximum load of 14MPa.As shown in Figure 1, the buoy was constructed from a cylin-drical pressure tank with an outer diameter measuring 76.2 cmand vertical length of 116.8 cm. The pressure vessel was later-ally sliced 23.5 cm from the top of the buoy and outfitted with abolted flange to provide access to the strain measurement systemon the interior of the buoy. The buoy was instrumented with 36strain gauges arranged in pairs that measured vertical and horizon-tal strain at 18 locations. The strain gauges were placed in threerings around the buoy and vertically spaced 16.5 cm apart, begin-ning with the top ring of 12 gauges located 34.7 cm from the buoytop (11.2 cm below the bottom of the buoy cap flange), a middlering located at 51.2 cm from the top (27.7 cm below the bottom ofthe buoy cap flange), and the bottom ring located 67.7 cm from thetop (44.2 cm below the bottom of the buoy cap flange). The gaugeswere uniformly distributed radially in 60 ° intervals (see Figure 1). Fig. 1 Dimensions of the buoy and strain gauge locations fromthe compression experiment. Strain gauges are placed in threerows (top, middle, and bottom) around the circumference ofthe buoy at ° intervals. Each point location has both a ver-tical and horizontal strain gauge. The colors correspond to thestrain results in Section 3. To simulate sea ice, the initial tank salinity was set to 27 ppt,which is slightly lower than typical Arctic sea surface salinity( ∼ ° C. The configuration and load condition resembles a uniaxialcompression test, illustrated in Figure 2. The sides of the ice sheetparallel to the loading direction were free of any confinement. Theside adjacent to the hydraulic press had a free slip boundary con-dition, but the opposite side was fixed. For the test, the hydraulicram load was gradually increased throughout the experiment untilreaching a peak value of approximately 14 MPa. The experimentlasted approximately 4 minutes with strain measurements recordedevery 0.5 seconds.
For each strain gauge location inFigure 1 in the horizontal (i.e., in the x-y plane tangent to thebuoy surface) and vertical strain was observed on the interior ofthe buoy. To relate this to the tractions acting on the exterior of thebuoy, we define a model predicting the buoy’s deformation givena particular traction field. This is done by solving the balance oflinear momentum equation with the assumption that inertial effectsare negligible. This is a reasonable assumption since the loadingrate of the ice on the buoy in the experiment is small. With thisassumption, the balance of linear momentum is given by ∇ · 𝜎 ( 𝑥 ) + 𝑏 ( 𝑥 ) = , (1)where 𝜎 ( 𝑥 ) is the Cauchy stress tensor and 𝑏 ( 𝑥 ) is the body forcevector on the buoy. Assuming a linear elastic constitutive model,2 nfrozen salt water H y d r au li c p i s t on Buoy Ice sheetDirection of loading 6.0 m . m Fig. 2 Schematic diagram of the setup for the buoy compres-sion experiment. The thick black lines indicate the test basinboundary and the red line indicates the boundary where the icewas frozen to the test basin wall. the stress-strain response follows Hooke’s law 𝜎 ( 𝑥 ) = 𝜇𝜀 𝑚𝑜𝑑 ( 𝑥 ) + 𝜆 tr ( 𝜀 𝑚𝑜𝑑 ( 𝑥 )) 𝐼, (2)where 𝜇 and 𝜆 are Lamé parameters and 𝜀 𝑚𝑜𝑑 is the model strain,which is defined as (assuming small strain) 𝜀 𝑚𝑜𝑑 ( 𝑥 ) = (cid:0) ∇ 𝑢 ( 𝑥 ) + (∇ 𝑢 ( 𝑥 )) (cid:124) (cid:1) , (3)where 𝑢 is displacement. We prescribe a Young’s modulus of 200GPa and a Poisson’s ratio of 0.3, which correspond to the materialproperties of the steel buoy used in the experiment. We then derivethe Galerkin weak form of Eq.(1) and arrive at the discretized finiteelement matrix equations 𝐾𝑑 = 𝑓 , (4)where 𝐾 is the stiffness matrix, 𝑑 is the nodal displacement vector,and 𝑓 is the external nodal force vector acting on the buoy surface Γ . For the purposes of inference, we only need to know predictedvalues of strain at the physical strain gauge locations on the buoy.Using the finite element basis functions that define the displace-ment 𝑢 with the definition of strain in Eq.(3), we are able to writethe strain at the observation locations as 𝜀 𝑚𝑜𝑑 = 𝐵𝐾 − 𝑓 , (5)where 𝐵 is a matrix that contains the derivatives of the finite el-ement shape functions. The size of the B matrix is constructedsuch that only the shape functions for the finite elements con-taining the location of the physical strain gauges are included. Itis a sparse matrix with 33 rows, one for each strain gauge, and203 ,
661 columns, one for each degree of freedom in the finite el-ement mesh. The nodal force vector 𝑓 can be further broken downinto 𝑓 = 𝐷 𝐴𝑝 (6)where 𝐷 is a matrix of shape function products integrated overthe external surface Γ , similar to a mass matrix, 𝑝 is a nodalpressure vector, and 𝐴 is a matrix that rotates and stamps thepressure vector from the buoy coordinate system, shown in Figure3, into the global coordinate system. The matrix 𝐴 has 203 , ,
073 columns, one for each degree of freedom we wishto infer. The pressure vector 𝑝 is composed of normal, 𝑝 𝑁 , andtangential components in the horizontal and vertical directions, 𝑝 𝐻 and 𝑝 𝑉 , respectively. The matrix 𝐴 transforms the results between Fig. 3 Schematic diagram of the unit vectors (horizontal: ^ t H ,vertical: ^ t V , and normal: ^ n ) used to project the pressure resultsbetween Cartesian coordinates and the buoy reference coordi-nates. The left image is a cross-sectional view looking from thetop of the buoy down, such that ^ t V is directed out of the pagetoward the reader. The right image is a cross-sectional view ofthe buoy rotated ° about the dashed line. (cid:174) r is the radius ofthe buoy. Cartesian global coordinates and the buoy coordinates through theexpression 𝑝 ( 𝑖 ) 𝑥 𝑝 ( 𝑖 ) 𝑦 𝑝 ( 𝑖 ) 𝑧 = (cid:104) 𝐴 ( 𝑖 ) 𝑛 𝐴 ( 𝑖 ) 𝐻 𝐴 ( 𝑖 ) 𝑉 (cid:105) 𝑝 ( 𝑖 ) 𝑁 𝑝 ( 𝑖 ) 𝐻 𝑝 ( 𝑖 ) 𝑉 , (7)where 𝐴 𝑛 = − ˆ 𝑛 ( 𝑖 ) − ˆ 𝑛 ( 𝑖 ) , 𝐴 𝐻 = ˆ 𝑡 ( 𝑖 ) 𝐻 ˆ 𝑡 ( 𝑖 ) 𝐻 , 𝐴 𝑉 = 𝑡 ( 𝑖 ) 𝑉 . (8)This linear elastic finite element model was implemented inFEniCS [33] using a tetrahedral mesh that was generated withthe CUBIT software, using a CAD geometry of the buoy used inthe experiment (Figure 1). The CAD files were provided by theUniversity of Washington Applied Physics Laboratory. We adopt a Bayesian approach to in-directly characterize the ice pressures 𝑝 𝑁 ( 𝑥 ) , 𝑝 𝐻 ( 𝑥 ) , 𝑝 𝑉 ( 𝑥 ) ap-plied to a defined region on the exterior surface of the buoy fromstrain gauge measurements on the interior of the buoy. Bayesianinference is a probabilistic framework that allows us to not onlyestimate the value of the pressures, but also to rigorously quantifyuncertainty stemming from noisy observations and the inabilityof a few observations to completely constrain the pressures overthe entire surface of the buoy. Let Γ 𝑝 ⊂ Γ denote the portion ofbuoy’s exterior surface that lies between 23.5 and 83.5 cm fromthe top of the buoy. Over this region, we want to infer the pressurefunction 𝑝 ( 𝑥 ) . We start by treating the external pressures appliedto the buoy and the observed strains as random variables, denotedby 𝑝 ( 𝑥 ) = [ 𝑝 𝑁 ( 𝑥 ) , 𝑝 𝐻 ( 𝑥 ) , 𝑝 𝑉 ( 𝑥 )] 𝑇 and 𝜀 𝑜𝑏𝑠 , respectively. Thecomponents 𝑝 𝑁 ( 𝑥 ) , 𝑝 𝐻 ( 𝑥 ) and 𝑝 𝑉 ( 𝑥 ) denote the stresses in thenormal, horizontal, and vertical directions, as defined in Figure 3.Our goal is to determine the conditional distribution of the pres-sures 𝑝 ( 𝑥 ) given the strain observations 𝜀 𝑜𝑏𝑠 . In general however, 𝑝 ( 𝑥 ) lies in an infinite dimensional function space, which makesinference more challenging. To overcome this, we instead look forpressure functions 𝑝 ( 𝑥 ) that lie within the span of a linear finiteelement basis; the same basis used above to represent the loadvector 𝑓 . This restriction allows us to characterize the pressurefunction 𝑝 ( 𝑥 ) with a finite number of coefficients and thus ap-ply standard Bayesian techniques. [34] provides more informationabout working directly in the infinite dimensional setting.3et 𝑝 (without the ( 𝑥 ) ) denote a vector containing nodal pres-sures. Our goal is then to characterize the posterior density 𝜋 ( 𝑝 | 𝜀 𝑜𝑏𝑠 ) . Bayes’ rule allows us to write this density as the prod-uct of a prior density 𝜋 ( 𝑝 ) and a likelihood function 𝜋 ( 𝜀 𝑜𝑏𝑠 | 𝑝 ) 𝜋 ( 𝑝 | 𝜀 𝑜𝑏𝑠 ) ∝ 𝜋 ( 𝜀 𝑜𝑏𝑠 | 𝑝 ) 𝜋 ( 𝑝 ) , (9)where the prior density models information known about the pres-sures before any observations are available, and the likelihoodfunction provides a way of comparing model predictions with theobservations. More information is provided on each of these com-ponents below. To obtain a prior distribution overthe finite dimensional vector 𝑝 , we first consider the continu-ous pressure function 𝑝 ( 𝑥 ) , which is spatially varying field. Itis natural to probabilistically describe the load function 𝑝 ( 𝑥 ) asa random field. In particular, we describe the prior distributionover the components 𝑝 𝑁 ( 𝑥 ) , 𝑝 𝐻 ( 𝑥 ) , and 𝑝 𝑉 ( 𝑥 ) with independentGaussian processes (see e.g., [35]). A Gaussian process definesa probability distribution over functions and can be interpreted asthe infinite-dimensional analog of the more common multivariateGaussian distribution. Like a multivariate Gaussian distribution,which is completely defined by its mean vector and covariancematrix, a Gaussian process is completely defined by a mean func-tion and a covariance kernel. The covariance kernel describes thecorrelation between loads at two points 𝑥 and 𝑥 (cid:48) while the meanfunction describes the average pressure at a single point. We use 𝜇 𝑁 ( 𝑥 ) , 𝜇 𝐻 ( 𝑥 ) , and 𝜇 𝑉 ( 𝑥 ) to denote the mean functions for eachcomponent and 𝑘 𝑁 ( 𝑥, 𝑥 (cid:48) ) , 𝑘 𝐻 ( 𝑥, 𝑥 (cid:48) ) , and 𝑘 𝑉 ( 𝑥, 𝑥 (cid:48) ) to denote thecorresponding covariance kernels. Together, these functions defineGaussian process distributions over each component 𝑝 𝑁 ( 𝑥 ) ∼ 𝐺𝑃 (cid:0) 𝜇 𝑁 ( 𝑥 ) , 𝑘 𝑁 ( 𝑥, 𝑥 (cid:48) ) (cid:1) (10) 𝑝 𝐻 ( 𝑥 ) ∼ 𝐺𝑃 (cid:0) 𝜇 𝐻 ( 𝑥 ) , 𝑘 𝐻 ( 𝑥, 𝑥 (cid:48) ) (cid:1) (11) 𝑝 𝑉 ( 𝑥 ) ∼ 𝐺𝑃 (cid:0) 𝜇 𝑉 ( 𝑥 ) , 𝑘 𝑉 ( 𝑥, 𝑥 (cid:48) ) (cid:1) . (12)We further assume that all of the covariance kernels take the form 𝑘 ( 𝑥, 𝑥 (cid:48) ) = 𝑘 (cid:32) tan − (cid:20) 𝑥 𝑥 (cid:21) , tan − (cid:34) 𝑥 (cid:48) 𝑥 (cid:48) (cid:35)(cid:33) 𝑘 (cid:16) 𝑥 , 𝑥 (cid:48) (cid:17) , (13)where 𝑥 𝑖 denotes component 𝑖 of the location 𝑥 , tan − (cid:104) 𝑥 𝑥 (cid:105) isthe angle around the buoy, 𝑘 is a 1d periodic covariance kerneldefining the meridional (horizontal) correlation of 𝑝 ( 𝑥 ) , and 𝑘 isa Matern kernel with 𝜈 = / √︁ 𝑘 ( 𝑥, 𝑥 ) ) is set to 4 . 𝑘 𝑁 and 0 . 𝑘 𝐻 and 𝑘 𝑉 . These standard deviations were chosento reflect the low probability that pressures would three times thesevalues due to the mechanical properties of the ice. The lengthscaleused in all meridional kernels is 𝜋 while the lengthscale used inall vertical kernels is 50 cm. These were chosen based on theanticipated smoothness of the pressure fields.As mentioned above, in order to use these Gaussian process de-scriptions with the finite element discretization described above,we need to discretize the pressure function 𝑝 ( 𝑥 ) . To do this withour Gaussian process prior distribution, we evaluate the mean func-tion at every FEM degree of freedom in Γ 𝑝 and evaluate the co-variance kernel for every pair of locations. The result is a finitedimensional Gaussian distribution that is more convenient to workwith than the infinite dimensional Gaussian process. The densityover the vector of nodal pressures 𝑝 is then given by 𝜋 ( 𝑝 ) = 𝑁 ( 𝜇 𝑝 , Σ 𝑝 ) = √︃(cid:12)(cid:12) 𝜋 Σ 𝑝 (cid:12)(cid:12) exp (cid:20) − (cid:0) 𝑝 − 𝜇 𝑝 (cid:1) 𝑇 Σ − 𝑝 (cid:0) 𝑝 − 𝜇 𝑝 (cid:1)(cid:21) , (14) where 𝑝 = (cid:34) 𝑝 𝑁 𝑝 𝐻 𝑝 𝑉 (cid:35) , 𝜇 𝑝 = (cid:34) 𝜇 𝑁 𝜇 𝐻 𝜇 𝑉 (cid:35) , Σ 𝑝 = (cid:34) Σ 𝑁 Σ 𝐻
00 0 Σ 𝑉 (cid:35) , (15)and 𝜇 𝑁 ,𝑖 = 𝜇 𝑁 ( 𝑥 ( 𝑖 ) ) , Σ 𝑁 ,𝑖 𝑗 = 𝑘 ( 𝑥 ( 𝑖 ) , 𝑥 ( 𝑗 ) ) , and 𝑥 ( 𝑖 ) denotes thelocation of mesh node 𝑖 in Ω 𝑝 . The Gaussian density 𝜋 ( 𝑝 ) denotesthe prior distribution over the external pressures. The likelihood function describesthe distribution of anticipated observations for a particular load 𝑝 ,i.e., What strains are likely to be observed for a particular load?Our finite element model is a relatively accurate representation ofthe buoy. We would therefore expect the modeled strain 𝜀 𝑚𝑜𝑑 ,given by Eq.(5)–(6), to be close to the observed strain 𝜀 𝑜𝑏𝑠 if thetrue pressures were used in the model. With an accurate model,the difference between 𝜀 𝑚𝑜𝑑 and 𝜀 𝑜𝑏𝑠 is then mostly a result ofnoise in the observations. We model this noise with a Gaussianrandom variable 𝑧 ∼ 𝑁 ( , 𝜎 𝑧 𝐼 ) and assume the additive form 𝜀 𝑜𝑏𝑠 = 𝜀 𝑚𝑜𝑑 + 𝑧. (16)The likelihood function 𝜋 ( 𝜀 𝑜𝑏𝑠 | 𝑝 ) is then a Gaussian density overthe difference 𝜀 𝑜𝑏𝑠 − 𝜀 𝑚𝑜𝑑 . In particular, 𝜋 ( 𝜀 𝑜𝑏𝑠 | 𝑝 ) = (cid:16) 𝜋𝜎 𝑧 (cid:17) − 𝑁 𝑧 / exp (cid:34) − 𝜎 𝑧 (cid:13)(cid:13)(cid:13) 𝜀 𝑜𝑏𝑠 − 𝐵𝐾 − 𝐷 𝐴𝑝 (cid:13)(cid:13)(cid:13) (cid:35) , (17)where the finite element model defined in Eq.(5) and (6) has beenintroduced to make the dependence on 𝑝 explicit. Multiplying the Gaussian prior inEq.(14) with the Gaussian likelihood in Eq.(17) yields the Bayesianposterior 𝜋 ( 𝑝 | 𝜀 𝑜𝑏𝑠 ) ∝ 𝜋 ( 𝜀 𝑜𝑏𝑠 | 𝑝 ) 𝜋 ( 𝑝 ) , which is also Gaussian anddefined by a mean vector 𝜇 𝑝 | 𝜀 and a covariance matrix Σ 𝑝 | 𝜀 . Tosimplify notation, consider the prior predictive covariance Σ 𝜀 = (cid:16) 𝐵𝐾 − 𝐷 𝐴 (cid:17) Σ 𝑝 (cid:16) 𝐵𝐾 − 𝐷 𝐴 (cid:17) 𝑇 + 𝜎 𝑧 𝐼, (18)and the Kalman gain 𝐺 = Σ 𝑝 (cid:16) 𝐵𝐾 − 𝐷 𝐴 (cid:17) 𝑇 Σ − 𝜀 . (19)The posterior mean and covariance are then given by 𝜇 𝑝 | 𝜀 = 𝜇 𝑝 + 𝐺 (cid:16) 𝜀 𝑜𝑏𝑠 − (cid:16) 𝐵𝐾 − 𝐷 𝐴 (cid:17) 𝑝 (cid:17) (20) Σ 𝑝 | 𝜀 = Σ 𝑝 − 𝐺 (cid:16) 𝐵𝐾 − 𝐷 𝐴 (cid:17) Σ 𝑝 . (21)The posterior mean 𝜇 𝑝 | 𝜀 describes the most likely pressures whilethe posterior covariance Σ 𝑝 | 𝜀 characterizes the remaining uncer-tainty. Note that the block structure in Σ 𝑝 and 𝐴 enable the pos-terior mean and covariance to be efficiently constructed withoutexplicitly building or storing the full covariance Σ 𝑝 . The Gaus-sianity of the prior and likelihood, combined with the use of a lin-ear model, allows us to compute the posterior analytically ratherthan having to use computationally expensive sampling methodslike Markov chain Monte Carlo.4 Results
We used a set of synthetic strain measurements that re-sulted from a known load configuration to verify that the inferenceframework works as expected and provides a reasonable character-ization of the applied pressures. The load configuration used forthe verification was composed of 3 rectangular patches that werelocated with respect to the 𝑥 -coordinate axis, which is aligned withthe 𝜃 = ◦ direction, as seen in Fig. 4. On the front sidewe applied a normal pressure of 4 MPa in a rectangular patch withdimensions 69.9 cm x 20.0 cm and on the back applied normalpressures of 2 MPa in two rectangular patches with dimensions34.4 cm x 20.0 cm, placed an equal distance on either side of the 𝑥 -axis. Fig. 4 shows the posterior mean of the normal load con-figurations. The inferred load pattern shows a region of elevatedpressure centered at 𝜃 = ◦ with a maximum pressure of 3.47 MPa.The pressures on the rear side of the buoy reach a maximum of1.65 MPa. The spatial extents of both front and rear regions ofhigh pressure compare well to the known applied load pattern. Thespatial pattern of the variances is indicative of the added informa-tion that the strain gauges provide and the differences in the lengthscales for the horizontal and vertical directions in the Gaussianprocess used to represent the inferred load fields. The variancetends to be high as one moves farther away from a strain gaugewith the decay in information being stronger in the vertical direc-tion. Interestingly, the prior lengthscale is actually shorter in thehorizontal direction. The slower increase of the posterior covari-ance in the horizontal direction is therefore an indication that thisstrain gauge configuration is more informative about horizontalvariations in the applied pressures.Figure 5 provides a more detailed look at the inferred load con-figuration for horizontal slices at vertical locations of 28.0 and35.25 cm below the buoy flange. There were aberrant negativeloads, but these negative excursions are very small and do not ex-ceed 0.87 MPa. The peak loads for the fore side of the buoy are3.31 MPa and located at 𝜃 = ◦ and on the aft side of the buoyare 1.56 MPa and located at 𝜃 = ◦ with the load configura-tion exhibiting a very slight asymmetry. There were small lobesat 𝜃 = ◦ and 260 ◦ , but these are relatively small in extent andmagnitude. As described inSection 2.2, we monitored the buoy’s strain using foil gauges thatwere affixed to the interior surface of the buoy as an ice sheetwas compressively loaded using a hydraulic ram. Figure 6, showshow both the vertically and horizontally aligned strain measuredat 18 locations on the buoy varied through the 300 second longexperiment. As the ice was initially loaded, the strains quicklyrise while the hydraulic ram applied approximately 7 MPa to theice until the 73 second mark. During this initial period, the ver-tical strains measured at the locations normal to the direction ofloading, 𝜃 = ◦ , were all tensile as one would expect fora hollow body pinched in the middle. In contrast, the horizontalstrains were mostly compressive.For time period between 73-251 seconds, the hydraulic ramwas set to its maximum output of 14 MPa through the duration ofthe experiment. During this timeframe, we observed a dramaticincrease in strain rate for all the locations, and most locationsreached a maximum value between 115 and 143 seconds. Therewere, however, a few exceptions to this strain trend. For example,the horizontal strains for the 𝜃 =
120 and 300 ◦ locations, whichface away and towards the hydraulic ram, respectively, continued toincrease after 115 seconds, but still exhibited a reduction in strainrate when the other locations reached their maximum strain values.In general, the strains were either purely tensile or compressivethrough the duration of the experiment. However, at the 𝜃 = ◦ locations, the vertical strains transitioned from compressionto tension, or vice versa, with the transitions occurring between135 and 161 seconds. In addition, the horizontal gauge at 𝜃 = ◦ transitioned from compression to tension at 124 seconds.The bottom horizontal strain gauge at 0 ° stopped collecting dataaround 10 seconds, therefore was removed from the dataset usedfor the load inversion. With the Bayesian inference frameworkintroduced in Section 2.4, we used the measured strains of the buoyinterior wall to infer the loads on the buoy’s exterior surface. Weperformed this inference on a subregion of the buoy surface, Γ 𝑝 ,that bounded the region that the ice was expected to be in contactwith the buoy. We independently processed each observation toobtain a time series of applied pressures, as shown in Figure 7and Figure 8. The inference results are presented as stress fieldsin the buoy coordinate system, with surface normal direction (pos-itive pointing towards the buoy interior), the horizontal direction(positive pointing clockwise) and the vertical direction (positivepointing up).For the initial 73 seconds, the stress levels on the buoy are verylow, which is to be expected since the strains were also small dur-ing this time interval. When the hydraulic ram output is increasedfrom 7 MPa to 14 MPa, the magnitude of the inferred stresses be-gin to increase until the horizontal, vertical, and normal stressesreach maxima of 2.7, 1.1, and 20.8 MPa, respectively, near a timeof 125 seconds into the experiment. The posterior mean stressfields at 100 and 125 seconds exhibit lobes near 0 ° that are con-sistent with the ice squeezing the buoy at this location and theformation of a crack at 0 ° . The horizontally aligned stress showsa line of zero stress at 𝜃 = ◦ with a positive and negative re-gion on either side of this zero stress line. This stress divergentpattern about the 0 ◦ line on the buoy is further evidence that acrack had formed at this location. In contrast, the vertical stresspattern shows a zero stress line aligned horizontally with a nega-tive stress (pushing down) above and a positive stress (pushing up)below this line, i.e., tangential stresses that are converging towardthis line. The largest vertical stresses are found along the bottomof the buoy’s flange, as shown at time 120s in Figure 9. Between100 and 125 seconds, the vertical stresses on the flange near 0 ° arepositive indicating that ice was pushing up on flange during thesetimes. As Figure 11 illustrates, this pressure on the flange agreeswith visible observations of the buoy being lifted by the flangeduring the experiment.Figure 7 provides a more detailed look at the inferred normalstresses along the middle ring of strain gauges (51.2 cm from thetop of the buoy) at 100, 120, 160, and 220 seconds. We see thelarge magnitude lobes pressure on either side of 0 ° and how theasymmetry increases with time favoring the sextant between 0 ° and60 ° . There are region of negative stress throughout the experimentthat we attribute to an unloading from a prestress in the buoy dueto the freezing in process, therefore we suspect that there is anunknown biases in the strain measurements that we were unable todetermine. The maximum inferred pressure in the lobe located inthe sextant between 0 ° and 60 ° was 20.5 MPa and the maximuminferred pressure in the lobe located in the sextant located between320 ° and 0 ° was 13.2 MPa . As asanity check, we calculate the posterior predictive mean strain atthe observation locations, i.e., the strain 𝜀 𝑝𝑜𝑠𝑡 , by using the pos-terior mean pressures as the boundary conditions for the FEM for-ward model. We compared these strain values against the actualobserved strain values, 𝜀 𝑜𝑏𝑠 , as a verification that the posteriorloads were matching the observations. The largest difference be-tween the mean posterior predictive and the observed strains was3 . × − , and occurred at 145 seconds, which corresponds toshortly after the time of highest observed strain.The inferred loads were also used to calculate the expected pos-terior displacement of the buoy surface. The largest displacementswere at 175 seconds, along the 0 ° axis (17 mm), where the buoydeformed inward. In addition to displacement in the direction of5 r on t B a ck Load Locations & Magnitudes Posterior Normal Loads Posterior Normal Variance
Fig. 4 Synthetic load configuration and inference results. ◦ ◦ ◦ ◦ ◦ ◦ -5.00.05.0MPa 28.00 cm Below Flange 0 ◦ ◦ ◦ ◦ ◦ ◦ -5.00.05.0MPa 36.25 cm Below Flange Fig. 5 Detail of inferred normal load configuration at 28.0 and36.25 cm below the flange. loading, the buoy bowed outward along the 90 ° and 270 ° axes(12 mm and 14 mm, respectively) in response to axial compres-sion. The displacement fields align well with the distribution ofthe normal loads, where highly compressive loads are found in thesame areas on the buoy where inward deformation occurred. Thedisplacement results at 200 seconds are shown in Figure 10 to il-lustrate the final deformation state of the buoy. These results helpto shape our understanding of the complex loading scenario thatcaused the observed strains. To evaluate howwell the inference framework worked for this experiment, we qual-itatively correlate the inferred load configurations with key eventsduring the physical experiment. Initially, the buoy was frozen intothe ice and had full contact between the buoy and the ice throughthe ice thickness. As the ice was compressively loaded, the buoyacted as an inclusion in the ice and induced localized ice crackingand deformation around the buoy. In some locations, the buoy de-coupled from the ice and transferred the load to highly localizedregions on the buoy. At approximately 100 seconds, a large crackformed in front of the buoy along the direction of loading alongwith two cracks off the backside of the buoy; one along the 180 ° axis and another emanating from the edge of the buoy near 240 ° .These two cracks on the backside separated a small segment of icefrom the rest of the ice sheet, as observed in images C and D ofFigure 11. The normal loads on the front of the buoy presented inFig. 8 showing compression in two lobes corroborates the pres-ence of a crack in the loading direction. The divergent horizontalloads in the 0 ° direction also indicate that the two sides of thiscrack were spreading apart as they pressed against the buoy. Thisis supported by observations that the crack became wider as theexperiment progressed.The ice did not remain planar throughout the experiment andstarted to buckle at 83 seconds, with some sections pushing upnearly 1 meter above the initial elevation. As it buckled, the ice6
100 200 300 − − ◦ − − ◦ − − ◦ − − ◦ − − ◦ Top HorizontalTop Vertical Middle HorizontalMiddle Vertical Bottom HorizontalBottom Vertical − − ◦ Time [s] S t r a i n [ m illi s t r a i n ] Fig. 6 Horizontal and vertical strain gauge measurements from the compression experiment through time. Each plot shows thetop, middle, and bottom gauges for one orientation on the buoy (dashed lines are vertical strain gauges, and solid lines are thehorizontal gauges). Note the scale difference for the ° and ° plots versus the others. ◦ ◦ ◦ ◦ ◦ ◦ -10.00.010.020.0MPa 0 ◦ ◦ ◦ ◦ ◦ ◦ -10.00.010.020.0MPa0 ◦ ◦ ◦ ◦ ◦ ◦ -10.00.010.020.0MPa 0 ◦ ◦ ◦ ◦ ◦ ◦ -10.00.010.020.0MPa
100 s 120 s160 s 220 s
Fig. 7 Rose plot showing the inferred normal load magnitudes(solid blue line) at times (a) 100, (b) 120, (c) 160, and (d) 220seconds. The standard deviations associated with these loadsis shown by the blue shaded region. The thick black line is 0.0MPa sheet lifted the buoy unevenly. The uneven lifting action was partlya result of the major crack that formed through the middle (at 0 ° )of the ice sheet. Wlipping between the ice and buoy was alsoobserved until the ice engaged with the buoy flange (see imageB in Fig 11). The inferred load configurations show large verticalloads on the flange and are highest upward load on the flange at 120seconds in Figures 8 and 9 consistent with the observations. Asdescribed previously, the ice load is higher in the sextant between0 ° and 60 ° matching the visual observations in the experimentwhere ice exhibiting more engagement with this sextant.Another key observation from the experiment was that the smallseparated section of ice on the back edge of the buoy rose with therest of the ice sheet as it buckled, which is an indicator that the asignificant portion of the buoy’s back side, i.e., between 90 ° and270 ° , was not in contact with ice for a large portion of the test(image C in Figure 11), which may explain why the normal loadswere inferred to be either small or even negative (see Fig. 7).Note that by inferring the loads acting on the buoy, we are notexplicitly modeling the ice-buoy contact forces. A loss of contactmanifests as a decrease in the load from the initial state of thebuoy (i.e., negative load in Fig. 7). The posterior probability ofa non-negative load could therefore be used to characterize theprobability of ice-buoy contact at any location on the buoy. Theseresults and corresponding observations illustrate how our inferenceframework could qualitatively capture the dynamic and complexloading conditions within the ice sheet that evolved during theexperiment. The good qualitative comparison between the ice behavior andinferred pressures shows promise for our approach. However, thereare certainly limitations as well. Certain aspects of the inferredload distribution, such as the high normal loads on the buoy inthe direction of loading, are easily explained with observationsfrom the uniaxial compression test. But, other features of theinferred loads are more difficult to justify using the visual obser-vations, such as the negative stress values on the buoy’s backside.Additional experiments with more careful control of the loadingconditions will be a better validation exercise for the approach presented in this paper. Another limitation is the small number ofobservation points ( 𝑁 𝑜𝑏𝑠 =
33) relative to the number of degreesof freedom ( 𝑁 𝑑𝑜 𝑓 𝑠 = , We have demonstrated a Bayesian inference framework that in-fers external load conditions on a buoy using a network of gaugesplaced on the interior surface of a steel cylindrical buoy. Using theinferred loads as inputs to a structural model of the buoy, we ver-ified the inferred results by comparing the modeled strain againstthe observed strains. The Bayesian inference approach presentedin this paper provides a robust framework for estimating the quan-tity of interest, in our case the most likely load configuration, butalso the associated uncertainty, which helps assess the quality ofthe inference results.Distributed observation of internal ice stress within the Arcticice pack is crucial for improving the understanding of sea ice dy-namics. To date, observations of sea ice internal stress are sparseor non-existent and certainly are not part of an ice monitoringsystem. We believe the framework described here could allow ex-isting buoy systems to be used as a network of stress sensors thatcan help capture the highly variable stress conditions in sea ice.Continuous, in situ internal stress measurements will allow us toobserve how these stresses develop, evolve, and propagate throughthe ice pack leading to improvements in our understanding of themechanics behind local processes like ridging and lead formation.This information will be crucial for improving sea ice models andproviding operational guidance in the region. The physical charac-ter and dynamic behavior of Arctic sea ice will continue to changein response to the warming climate and instrumentation that canprovide new insight into the mechanical behavior sea ice will beessential. We believe the Bayesian inference framework presentedin 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 PhysicsLaboratory for providing the buoy used in the experiment andthe strain gauge data. We would also like to thank Dr. AndrewDavis for taking the time to provide constructive feedback for thismanuscript.8 his study was supported by the Office of Naval ResearchArctic and Global Prediction program (code 32) under awardN0001418MP00501.
References [1] Comiso, J. C., 2012, “Large Decadal Decline of the Arctic Multiyear Ice Cover,”Journal of Climate, (4), pp. 1176–1193.[2] Screen, J. A. and Simmonds, I., 2010, “The Central Role of Diminishing SeaIce in Recent Arctic Temperature Amplification,” Nature, (7293), p. 1334.[3] Jeffries, M. O., Overland, J. E., and Perovich, D. K., 2013, “The Arctic,” Phys.Today, (10), p. 35.[4] Carmack, E., Polyakov, I., Padman, L., Fer, I., Hunke, E., Hutchings, J., Jackson,J., Kelley, D., Kwok, R., Layton, C., et al., 2015, “Toward Quantifying theIncreasing Role of Oceanic Heat in Sea Ice Loss in the New Arctic,” Bulletinof the American Meteorological Society, (12), pp. 2079–2105.[5] Strong, C. and Rigor, I. G., 2013, “Arctic Marginal Ice Zone Trending Widerin Summer and Narrower in Winter,” Geophysical Research Letters, (18), pp.4864–4868.[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, , pp. 300–317.[7] Ho, J., 2010, “The Implications of Arctic Sea Ice Decline on Shipping,” MarinePolicy, (3), pp. 713–715.[8] Melia, N., Haines, K., and Hawkins, E., 2016, “Sea Ice Decline and 21st Cen-tury Trans-Arctic Shipping Routes,” Geophysical Research Letters, (18), pp.9720–9728.[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 CanadianArctic,” Geophysical Research Letters, (23).[10] Stroeve, J., Markus, T., Boisvert, L., Miller, J., and Barrett, A., 2014, “Changesin Arctic Melt Season and Implications for Sea Ice Loss,” Geophysical ResearchLetters, (4), pp. 1216–1225.[11] Schreyer, H., Sulsky, D., Munday, L., Coon, M., and Kwok, R., 2006, “Elastic-Decohesive Constitutive Model for Sea Ice,” Journal of Geophysical Research:Oceans, (C11).[12] Hutchings, J. K., Roberts, A., Geiger, C. A., and Richter-Menge, J., 2011,“Spatial and Temporal Characterization of Sea-Ice Deformation,” Annals ofGlaciology, (57), pp. 360–368.[13] Marsan, D., Stern, H., Lindsay, R., and Weiss, J., 2004, “Scale Dependenceand Localization of the Deformation of Arctic Sea Ice,” Physical review letters, (17), p. 178501.[14] Rampal, P., Weiss, J., Marsan, D., Lindsay, R., and Stern, H., 2008, “ScalingProperties of Sea Ice Deformation from Buoy Dispersion Analysis,” Journal ofGeophysical Research: Oceans, (C3).[15] Johnson, J., Cox, G., and Tucker, W., 1985, “Kadluk Ice Stress MeasurementProgram,” Proceedings of the Eighth International Conference on Port andOcean 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 DriftingPack Ice,” Cold regions science and technology, (2), pp. 119–139.[18] Coon, M. D., Knoke, G. S., Echert, D. C., and Stern, H. L., 1993, “Contem-poraneous Field Measurements of Pack Ice Stress and Ice Strain Measurementsfrom SAR Imagery,” OCEANS’93. Engineering in Harmony with Ocean. Pro-ceedings , IEEE, pp. III31–III36.[19] Richter-Menge, J. and Elder, B., 1998, “Characteristics of Pack Ice Stress in theAlaskan Beaufort Sea,” Journal of Geophysical Research: Oceans, (C10),pp. 21817–21829.[20] Hata, Y. and Tremblay, L., 2015, “Anisotropic Internal Thermal Stress in SeaIce From the Canadian Arctic Archipelago,” Journal of Geophysical Research:Oceans, (8), pp. 5457–5472.[21] Zweng, M. M., Boyer, T. P., Baranova, O. K., Reagan, J. R., Seidov, D., andSmolyar, I. V., 2018, “An Inventory of Arctic Ocean Data in the World OceanDatabase,” Earth System Science Data, (1), p. 677.[22] Cox, G. F. and Johnson, J. B., 1983, “Stress Measurements in Ice,” USA ColdRegions Research and Engineering Laboratory, Hanover, NH.[23] Metge, M., Strilchuk, A., and Trofimenkoff, P., 1975, “On Recording Stresses inIce.” Proc. IAHR Third International Symposium on Ice Problems , USA ColdRegions Research and Engineering Laboratory, Hanover, New Hampshire, pp.459–468.[24] Templeton, J., 1979, “On Measurement of Sea Ice Pressures,”
Proc. Fifth Inter-national Conference on Port and Ocean Engineering under Arctic Conditions ,Norwegian Institute of Technology, Trondheim, Norway, pp. 73–88.[25] Johnson, J. and Cox, G., 1980, “The OSI Ice Stress Sensor,”
Proceedings ofthe Workshop on Sea Ice Field Measurements , Vol. 80, Center for Cold OceanResources Engineering St. John’s, pp. 193–207.[26] GeoKon, 2018,
Instruction Manual: Model 4350BX Biaxial Stressmeter ,GeoKon, Lebanon, New Hampshire.[27] Borcea, L., 2002, “Electrical Impedance Tomography,” Inverse problems, (6),p. R99.[28] Schulz, R. B., Ale, A., Sarantopoulos, A., Freyer, M., Soehngen, E., Zien-tkowska, M., and Ntziachristos, V., 2010, “Hybrid System for SimultaneousFluorescence and X-Ray Computed Tomography,” IEEE transactions on medi-cal imaging, (2), pp. 465–473. [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 froma Full-Stokes and Simpler Models for Pine Island Glacier, West Antarctica,”Geophysical Research Letters, (14).[30] Steiner, B., Saenger, E. H., and Schmalholz, S. M., 2008, “Time Reverse Mod-eling of Low-Frequency Microtremors: Application to Hydrocarbon ReservoirLocalization,” Geophysical Research Letters, (3).[31] Vanik, M. W., Beck, J. L., and Au, S., 2000, “Bayesian Probabilistic Approachto Structural Health Monitoring,” Journal of Engineering Mechanics, (7),pp. 738–745.[32] Brucker, L., DIinnat, E., and Koenig, L., 2015, “Aquarius L3 Weekly Polar-Gridded Brightness Temperature and Sea Surface Salinity, Version 5,” USANASA 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 FEniCSProject Version 1.5,” Archive of Numerical Software, (100).[34] Stuart, A. M., 2010, “Inverse Problems: a Bayesian Perspective,” Acta Numer-ica, , pp. 451–559.[35] Rasmussen, C. E. and Williams, C. K. I., 2005, Gaussian Processes for MachineLearning , The MIT Press.[36] Cui, T., Law, K. J., and Marzouk, Y. M., 2016, “Dimension-independentlikelihood-informed MCMC,” Journal of Computational Physics, , pp. 109–137.[37] Cotter, S. L., Roberts, G. O., Stuart, A. M., and White, D., 2013, “MCMCmethods for functions: modifying old algorithms to make them faster,” Statisti-cal Science, pp. 424–446. ig. 8 Inferred loads across the buoy looking down the direction of loading ( ° ) at times 100, 120, 160, and 220 seconds, respec-tively. The first column shows the horizontal loads, the middle column shows the normal loads, and the last column the verticalloads. The color bar units are MPa. ig. 9 Detailed view of vertical pressure on the flange at 120slooking from below. Positive quantities are pointing into thepage. There is a clear lift corresponding to the visible flange-ice contact in image B of Figure 11. Units are in MPa. o o Fig. 10 Calculated displacements at 200 seconds using the in-ferred load results. The buoy is warped by a scaling factor of2.0, and is show oriented along the ° and ° axes ) 45 o T=0 B) 45 o T=125C) 225 o T=150 D) 160 o T=250