1. Introduction
Accelerated melting of glaciers and ice caps has raised serious concerns about sea level rise. As we work towards a solution to address these concerns, it has become a chief priority to rapidly improve predictions of future changes in global ice mass balance. In particular, numerical simulations projecting ice loss have uncovered a strong sensitivity to mechanical and/or rheological weakening of the shear margins of streaming ice (Joughin and others, Reference Joughin2014). To accurately project sea level rise, future models will require careful treatment of shear margins. This necessitates a deeper understanding of the flow dynamics at shear margins and how streaming flow relates to the constitutive flow law for ice (Glen, Reference Glen1955, Reference Glen1958).
Streaming ice is characterized by strong shearing, wet-based beds with significant basal sliding and/or deformable till and the development of anisotropic crystal fabric (Jackson and Kamb, Reference Jackson and Kamb1997). Streaming flow occurs in two major groups of glaciers: polythermal and temperate. Regions of streaming flow are responsible for draining the major ice sheets (Azuma, Reference Azuma1994) and alpine regions (Lliboutry and Duval, Reference Lliboutry and Duval1985). There remains a lack of systematic in-situ studies in natural settings of streaming flow, and most in-situ (e.g. Faria and others, Reference Faria2014) and in-lab (Nye, Reference Nye1953; Glen, Reference Glen1955) studies of natural ice have focused on bodies of ice with frozen beds that experience minimal shear strain across planes normal to both the bed and margins, and where the viscosity at the base of the ice controls most of the velocity field.
In an effort to bridge this knowledge gap, we seek to contribute to kinematic measurements of streaming ice through instrumentation. Our instrument is a sensor, containing a 3-axis accelerometer and a 3-axis magnetometer, that measures tilt. We drilled multiple boreholes close to the shear margin of Jarvis Glacier in Alaska and installed our tilt sensor systems in two boreholes, collecting around a year of uninterrupted tilt data from Aug. 2017–Aug. 2018. From the tilt sensors' gravimetric and magnetic data, we calculated borehole deformation by tracking the orientation of the sensors over time. We also evaluated our results against flow dynamic parameters derived from Glen's flow law and explored the parameter space of the stress exponent n and enhancement factor E. When tuning n, we referenced the experimental n values of Goldsby and Kohlstedt (Reference Goldsby and Kohlstedt1997a,Reference Goldsby and Kohlstedtb, Reference Goldsby and Kohlstedt2001) which they have correlated with stages of creep. When tuning E alongside n, we calculated the mismatch value between theoretical models (sets of n and E) and our observed results to generate a mismatch index over a range of n and E.
2. Methods
2.1. Tilt sensor
To obtain kinematic measurements of streaming ice, we developed an inexpensive tilt sensor built from accessible parts that can be installed at multiple depths in a borehole. A simple sensor schematic is shown in Fig. 1 and highlights the core components. Prior studies using borehole inclinometers guided our sensor development process, particularly on optimizing sensor length to prevent strain analysis complications (Hooke and Hanson, Reference Hooke and Hanson1986) and maintaining a large enough (>4) aspect ratio for tracking the sensor. For an inclinometer embedded in a Newtonian or Glen-type fluid to behave like a Lagrangian material line element, a large enough aspect ratio was needed to ensure the inclinometer was primarily impacted by horizontal shear strain instead of vertical strain when embedded in ice (Jaber and others, Reference Jaber2013).
2.2. Tilt sensor system
The design concept for our tilt sensor system is a string of connected tilt sensors, each weighting less than a pound, communicating along a common data bus. The sensor string is lowered into the borehole and frozen-in to collect gravity and magnetic data, which are transmitted up to a Campbell Scientific CR-1000 data logger at the surface. The system is powered by two 12 V sealed lead acid batteries and all system components are connected by 22-gauge 4-conductor shielded cables (with up to 10 m spacing between adjacent components) containing a ground wire for redirecting short circuits. To eliminate significant stress on the connection between the cable and the connector, we implemented a strain relief system by making cable loops from both ends of each tilt sensor and hose clamping the loops to the sensors to redirect the load onto the loops.
Communication and data transmission for our tilt sensor system uses two serial communication standards, transistor-transistor logic (TTL) and RS-485. The LSM303C breakout board and Arduino Pro Mini within each sensor communicate using TTL, while longer-distance communication (up to 80 m) between the sensors and data logger at the surface necessitates the use of RS-485. We installed a linear LTC1484 RS-485 transceiver in each tilt sensor and an SDM-SIO1A I/O expansion module in the data logger to convert between TTL and RS-485 signals.
Data collection for our tilt sensor system uses an addressable bus approach. Each tilt sensor is coded with a unique address number and programmed to respond when its address number is called. At 6-h intervals, the data logger powers on and sends out a series of call commands to the sensors to trigger transmission of the sensors' current data over the data bus to the data logger. The call commands are sent out staggered in time to ensure that only one Arduino Pro Mini is active on the communication bus at any time to prevent data congestion. We chose a 6-h sampling time for each sensor as a trade-off between dense temporal sampling (to minimize integration error in our sensor fusion algorithms) and limited data storage on the data logger. Power supply problems (Keller and Blatter, Reference Keller and Blatter2012) turned out to not be an issue, as after 18 months of data recording the batteries were down to ~11.2 V from 13.2 V fully charged. Our sampling rate was dictated by storage space in the data logger rather than power consumption. As the data size transmitted through the system at any one time was small, we were able to use a slower 9600 baud rate for the Arduino Pro Mini and data logger, which helped protect against communication errors over the long cable.
2.3. Sensor calibration
We performed on-site magnetic calibration of our sensors prior to installation in the boreholes on Jarvis, to correct for magnetometer interference. We positioned each sensor at an inclination of 0° (upright) in a calibration jig and took eight sets of magnetic x and z readings, two sets every 90° as we rotated the sensor 360° about its longitudinal axis. We then plotted each sensor's readings on a scatter chart (representing the xz-plane) and fitted the readings with a direct least-squares circle, allowing us to determine the x and z offsets of the circle from the origin (0, 0) and apply the magnetic corrections. The sensors appear well calibrated, as our observed 0 G offsets (− 0.135 G averaged over all sensors' x and z readings, deviating by around ± 0.05 G when separating out offsets between boreholes and x/z readings) are well within the LSM303C datasheet's (available athttps://www.st.com/resource/en/datasheet/dm00089896.pdf) quoted typical0 Goffset accuracy of ±1 G.
We did not perform gravity calibrations as it has no effect on our analysis since we only think about relative orientation and are calculating changes in angles rather than absolute angles. Any 0-G offset for each sensor would apply throughout its recording and would not affect our calculations of the net angle of rotation ψ of the sensor (detailed in section 2.5). We did however perform lab tests for our sensors over a range of ψ [0° 180°] and inclinations θ [0° 360°] to ensure that our sensor outputs matched expected values. These lab tests are not detailed gravity calibrations but rather checks as they were performed without a real calibration jig to get the sensors oriented in precisely the expected angles.
2.4. Study area
Jarvis Glacier is a polythermal mountain glacier (Miner and others, Reference Miner2018) in the Eastern Alaskan Range, which has streaming ice and a simple geometry. Jarvis is 5.2 km2 in surface area, 8 km long and up to 220 m thick, which is sufficiently shallow to allow drilling of multiple surface-to-bed cores in a single season. It has a headwall elevation of 2000 m and a small tributary glacier joins the main glacier from the southeast. The air temperatures on Jarvis range from −20°C to 10°C and relative humidity ranges from ${25}\percnt$ to ${90}\percnt$. Satellite imagery of Jarvis is shown in Fig. 2.
We drilled multiple boreholes at our research site on Jarvis using an Ice Drilling Design and Operations Badger-Eclipse electromechanical drill and installed our tilt sensor systems in three boreholes. The first borehole lies closest to the glacier centerline at site JA (63.4750° N, −145.6753° W, 1621 m) and is 80 m deep. The second borehole is at the shear margin close to the valley walls at site JE (63.4743° N, −145.6766° W, 1625 m) and is 18 m deep. We had installed our tilt sensor system in a third borehole at site JC down-glacier of JE, but lost our instruments due to crevasse formation. We spaced the sensors more closely with depth in the borehole to capture in higher resolution the greater rate of change of borehole tilt with depth, since basal resistance and overburden pressure causes deformation to be greatest near the bed. We installed 13 sensors in JA at 2.00, 10.35, 19.65, 28.40, 37.15, 45.95, 54.50, 58.70, 62.85, 67.00, 69.15, 70.35 and 71.70 m and seven sensors in JE at 3.50, 7.58, 10.84, 13.09, 15.39, 16.63 and 17.84 m.
2.5. Calculating borehole deformation and flow dynamics
From the tilt sensors' gravimetric and magnetic data we calculated the sensor's tilt inclination θ and azimuth ϕ to inform us of the sensor's orientation (Appendix A), and recreated the orientation of the sensors in a Cartesian coordinate system using three-dimensional rotation. We then tracked the orientation of the sensors and their rotation axes over time, and calculated borehole deformation using the net angle of rotation ψ of the sensor in the down-glacier direction within a two-dimensional down-flow xz-plane as illustrated in Fig. 3. From ψ we calculated the flow dynamic parameter shear strain rate $\dot {\varepsilon }$ (Appendix B) and vertical gradients of down-glacier velocity du/dz. Since the shear strain rate is defined as one half the tangent of the rate of decrease in an initially right angle per unit time, we multiplied $\dot {\varepsilon }$ by two to obtain du/dz. Our two-dimensional analysis is restricted to du/dz measurements in the down-flow xz-plane. We lacked control on du/dy measurements as our sensors were not installed at similar depths between boreholes, and we had negligible du/dx measurements as our sensors remained nearly vertical through the data collection period, resulting in not much Δx to measure over.
3 Results
3.1. Data description
We collected around a year of uninterrupted tilt data from Aug. 2017–Aug. 2018 from the boreholes at sites JA and JE. Excluding sensors that had melted out onto the surface and/or were not coupled to the borehole wall, we tracked the orientation of the remaining sensors and their rotation axes over time. The rotation axis for each time step was calculated from the cross product of a sensor's orientation at an initial time t1 and at a subsequent time t2 following a time step. Polar plots tracking the sensors' rotation axes over time are shown in Fig. 4. The rotation axes of JA sensors travel back and forth systematically along arced paths that trend towards and cluster around 215° azimuth, indicating a general down-glacier flow (the down-flow azimuth relative to magnetic north is ~305° and by definition the corresponding down-flow rotation axis is 90° counterclockwise of the down-flow azimuth). The rotation axes of JE sensors travel along similar systematic arcs that trend past 215° without any localized clustering, which indicates an element of a cross-glacier flow along with the general down-glacier flow.
Our projection of the sensors onto the down-flow xz-plane to calculate ψ assumes a general down-glacier flow and that has held true for JA and JE, with an element of the cross-glacier flow in JE which might result in underestimation of its two-dimensional kinematic measurements and subsequent derived calculations of flow dynamics. As streaming ice is characterized by significant basal sliding, our calculations of an internal deformation flow from the tilt data constitute only a small fraction of the total velocity (Harper and others, Reference Harper2001). Factoring in a basal sliding component, the cross-glacier flow observed in JE as well as the path of the sensors' rotation axes (some reversing over time, though most trended towards 215° by the end of the data collection period) could be necessary adjustments for sliding over a rough bed and/or seasonal adjustment from surface or basal mass loss. Seasonal variations in the basal meltwater input and drainage system morphology/efficiency can change the stress partitioning between the bed and lateral margins; we plotted a seasonal time series of the net angle of rotation ψ of the sensors in Fig. 5 but observed no discernible seasonal variations in the rate of tilt. We also considered the effects of vertical shortening and longitudinal strain on our measurements, but found them negligible due to the small tilt inclination (often significantly less than 10°) of the sensors that persisted through the data collection period. We calculated borehole deformation and vertical gradients of down-glacier velocity du/dz as detailed in section 2.5 and show our du/dz results in Fig. 6.
3.2. Error analysis
We resolved the two sources of uncertainty, instrumental and field, for our tilt sensors in this study to generate error bars with ${95}\percnt$ confidence intervals (two standard deviations) for our results.
Our tilt sensor has an inherent instrumental uncertainty determined by its precision, based off the limitations of the analog-to-digital converter used by the sensor. The LSM303C has multiple linear acceleration full scales (FS), and we used FS = ±2 G since we are detecting small accelerations. The LSM303C datasheet indicates a sensitivity of 0.061 mG per least significant bit (LSB) at FS = ±2 G for the accelerometers and 0.58 mG/LSB at FS= ±16 G for the magnetometers. For temperature dependence, the datasheet indicates a linear acceleration sensitivity change vs temperature value of ${0.01}\percnt /^{\circ }{\rm C}$, which is small enough to be negligible. The datasheet also indicates that root-mean-square (RMS) noise is 1 mG for the accelerometers and 3.5 mG for the magnetometers. Based off our lab data, the linear acceleration observed noise is ~0.7 mG and the magnetic observed noise is ~3.2 mG, both values close to and smaller than their corresponding datasheet numbers. To calculate instrumental uncertainty, we collected tilt data from two-parallel clamped sensors recording for 7 days in a regulated constant lab setting and ran the data through the algorithms detailed in section 2.5 to calculate inclination and azimuth, taking the standard deviation for each parameter.
Each tilt sensor also has a unique field uncertainty, determined by the borehole location and installation depth. As our sensors are connected by cables, we also considered the effect of top sensors melting out onto the surface on the accuracy of measurements of the remaining sensors. We observed the tilt of the remaining sensors in the same time frame that the top sensors melted out but noted no significant fluctuations of their measurements. We also considered temperature effects particularly for sensors installed not far below the glacier surface, but all sensors above 10 m had melted out and we found no discernible correlation when comparing the remaining sensors' field uncertainties with their corresponding seasonal temperature fluctuations. To calculate field uncertainty for each sensor, we determined their total uncertainty for inclination and azimuth from the standard deviation of their residuals derived from our observed and fitted (by running our observed data through a Butterworth low-pass filter) data and subtracted the corresponding instrumental uncertainty.
As the flow dynamic parameters are derived from inclination and azimuth, we performed error propagation to determine the derived uncertainties up to du/dz. We performed a Monte Carlo simulation by generating 1000 normally distributed iterations of inclination and azimuth sets per sensor per time step and ran each set through the algorithms detailed in section 2.5, using the standard deviations for the iterated flow dynamic parameters to generate the error bars shown in Figs. 6 and 7. The observed du/dz instrumental uncertainty is ~7.50 × 10−4 a−1 (close to and smaller than the RMS instrumental uncertainty of ~0.001 a−1) and the observed du/dz field uncertainties are highly variable between sensors due to the different environmental conditions as detailed above, with an average value of ~0.0085 a−1.
4. Discussion
4.1. Tuning the stress exponent n
We plotted the du/dz results at JA and JE against theoretical solutions to Glen's Flow Law ina laminar flow (Cuffey and Paterson, Reference Cuffey and Paterson2010) for a range of n values in Fig. 7, using the physical parameters measured in the borehole and setting the enhancement factor E to 1 for isotropic ice since we are unable to account for and separate the collective effects of variables such as fabric and water content that all lead to deviations from unity.
The results at JA fit closest to n = 2.9, which is close to the empirical fit value of 3 for glacier ice and from the ice deformation experiments of Glen (Reference Glen1955). The results at JE fit closest to n = 3.4, which falls short of the value of ~4 found for dislocation creep from the ice deformation experiments of Goldsby and Kohlstedt (Reference Goldsby and Kohlstedt1997a,Reference Goldsby and Kohlstedtb, Reference Goldsby and Kohlstedt2001). Our range of best fit n for JA and JE suggests that the creep regime of the ice on Jarvis lies between basal slip (n = 2.4) and dislocation creep (n = 4) Goldsby and Kohlstedt (Reference Goldsby and Kohlstedt1997a,Reference Goldsby and Kohlstedtb, Reference Goldsby and Kohlstedt2001). We also note that the best fit n values for JA and JE are higher than expected in regions of low stress (defined as σ ≤0.1 MPa, max. σJA ~78 kPa and max. σJE ~18 kPa), as Goldsby and Kohlstedt (Reference Goldsby and Kohlstedt1997a,Reference Goldsby and Kohlstedtb, Reference Goldsby and Kohlstedt2001) had experimentally determined that n approaches 2.4 in regions of low stress. The higher n values could be attributed to the observed high total strains (max. $\varepsilon ^{t}_{\rm {JA}}\; {\sim }{410}\percnt$ and max. $\varepsilon ^{t}_{\rm {JE}}\; {\sim }{850}\percnt$): The high strain rates coupled with potential dynamic recrystallization cause anisotropic development as the ice crystal fabric is increasingly aligned along preferred orientations, creating planes of weakness that speed up ice flow (Van der Veen and Whillans, Reference Van der Veen and Whillans1990; Wilson and Peternell, Reference Wilson and Peternell2011). We however acknowledge our incomplete understanding of Jarvis' strain history, so the period of total strain build-up at the boreholes is unknown.
4.2. Exploring the parameter space of n and E
It is reasonable to assume that there is anisotropic development of Jarvis ice at the shear margins causing the enhancement factor E to deviate from the isotropic value of 1 since for fabrics, E depends on the orientation of the deformation. We observed in Fig. 7 that tuning n affects the curvature of our theoretical solutions, though their fit with our observed results could be improved by simultaneously tuning E to shift the theoretical models laterally. We tuned both n and E (with a theoretical upper E limit of 10 for anisotropic ice with single-maximum fabric) in our theoretical models (sets of n and E) and calculated their mismatch values relative to our observed results using the root sum square. The du/dz mismatch indexes for JA and JE are shown in Fig. 8. We acknowledge that the assignment of an enhancement factor to the constitutive law is an expedient technique to quantify the variations of the strain rate not explained by stress and temperature, and so have avoided attempts to parse the physical meaning of E since we are unable to account for and separate the collective effects of variables such as fabric and water content that all lead to deviations from unity.
The JA mismatch index has a best fit at n ≈ 2.9 when E = 1, which decreases to n ≈ 2.7 as E increases to 10. The JE mismatch index has a best fit at n ≈ 3.4 when E = 1, which decreases to n ≈ 3.2 as E increases to 10. We note that n has a greater impact on the fit than E, which is not surprising given the difference between an exponent and scaling factor. In our discussion we have assumed a degree of homogeneity of the borehole physical parameters with depth to keep n and E constant throughout the borehole, as tuning n and E for different depths makes interpretation increasingly complex. Keeping n and E constant throughout the borehole allows us to describe the ice with simple parameters, though we acknowledge that simplifying n and E as local factors makes them valid only for the specific set of stress-temperature-fabric conditions in the borehole; n and E can and likely do change all the way down the borehole.
5. Conclusions
We have developed an inexpensive tilt sensor built from accessible parts that has delivered significant cost savings and is a viable substitute to commercial sensors at ~20% the cost. To test the viability of our sensors we installed our tilt sensor systems in two boreholes drilled close to the shear margin of Jarvis Glacier, Alaska to obtain kinematic measurements of streaming ice, and successfully collected around a year of uninterrupted tilt data. The development of our tilt sensor system and associated processing algorithms is open source, improving the accessibility of borehole geophysics.
From the tilt sensors' gravimetric and magnetic data we calculated borehole deformation by tracking the orientation of the sensors over time, which brought out interesting kinematic differences between the boreholes at sites JA and JE. The flow in JA, which lies closer to the centerline, displays a general down-glacier trend while the flow in JE, which is at the shear margin, displays an element of the cross-glacier flow with a general down-glacier trend which might result in underestimation of its two-dimensional kinematic measurements and subsequent derived calculations of flow dynamics. As streaming ice is characterized by significant basal sliding, our calculations of the internal deformation flow from the tilt data constitute only a small fraction of the total velocity; the cross-glacier flow observed in JE could be necessary adjustments for sliding over a rough bed and/or seasonal adjustment from surface or basal mass loss.
Beyond kinematic measurements, our brief investigation into the flow dynamics of streaming ice on Jarvis yielded higher values of the stress exponent n, particularly at the shear margin, than are expected in regions of low stress. The higher n values could be attributed to the observed high total strains coupled with potential dynamic recrystallization, causing anisotropic development and consequently sped up ice flow. Jarvis' n values place the creep regime of the ice between basal slip and dislocation creep Goldsby and Kohlstedt (Reference Goldsby and Kohlstedt1997a,Reference Goldsby and Kohlstedtb, Reference Goldsby and Kohlstedt2001). We were also able to explore the parameter space of n and the enhancement factor E to give us with a range of best fit parameters for our solution. A constant n and E throughout the borehole allows us to describe the ice with simple parameters, though we acknowledge the limitations of simplifying n and E as local factors (n and E can and likely do change all the way down the borehole). We found that tuning E towards a theoretical upper limit of 10 for anisotropic ice with single-maximum fabric reduces the n values by 0.2.
Data
The data reported in our manuscript are available in the Arctic Data Center at DOI: 10.18739/A2348GG12 or on request.
Acknowledgments
We thank several people who were essential to this project: D. Collins, C. Deck, M. Waszkiewicz and co-PIs P. Koons and K. Kreutz. This project was funded by the US National Science Foundation grant OPP 1503653. Logistical support was provided by CH2M HILL Polar Services, Talkeetna Taxi and 40-Mile Air.
Appendix A. Inclination and Azimuth
We first transform our gravimetric and magnetic data, which are in left-handed Cartesian coordinates, to more familiar right-handed Cartesian coordinates. The transformation is shown in Fig. 9 and produces the right-handed gravity vector A and magnetic vector M.
Inclination θ is the angle from the vertical gravity axis A from which the sensor's vertical axis Z′ has tilted, as illustrated in Fig. 10. We can find θ by calculating the angle between A and Z′ using the dot product
Azimuth ϕ is the direction of the sensor's vertical axis Z′ relative to magnetic north (0° azimuth). To measure the azimuth on a plane normal to the vertical gravity axis A, we project Z′ and the magnetic vector M onto the plane normal to A as illustrated in Fig. 11. The Z′ projection is Z′⊥ A
and the M projection is M′, representing the component of the magnetic north vector on the plane normal to A
Z′⊥ A and M′ lie in the same plane and are normal to A. We can find ϕ by calculating the angle of rotation from Z′⊥ A to M′ about A using the dot product
To determine the direction of rotation of ϕ, we calculate the dot product between (Z′⊥ A) × M′ and A. The cross product (Z′⊥ A) × M′ yields a vector that will be parallel, anti-parallel or perpendicular to A. A positive dot product indicates (Z′⊥ A) × M′ and A are parallel and ϕ is being calculated in our intended counterclockwise rotation direction. A negative dot product indicates (Z′⊥ A) × M′ and A are anti-parallel and ϕ is being calculated in the clockwise rotation direction, so we subtract 360° from ϕ to flip the rotation direction to counterclockwise. A dot product of zero indicates (Z′⊥ A) × M′ and A are perpendicular, meaning Z′⊥ A and M′ lie on top of each other on the plane normal to A which makes ϕ = 0°.
Appendix B. Shear Strain Rate
Glen's Flow Law (Glen, Reference Glen1955) defines the relation of strain rate $\dot {\varepsilon }$ to deviatoric stress σ. In our two-dimensional analysis, we restrict measurements of the sensor's tilt to the down-flow xz-plane and are only interested in the tensor component $\dot {\varepsilon }_{zx}$
where E is the enhancement factor, A is the rate factor in an isothermal case determined by the average borehole temperature over the data collection period and n is the stress exponent. The down-flow shear stress component is σzx = ρgzsinα, where z is borehole depth, ρ is density of glacier ice (917 kg/m3), g is gravitational constant (9.81 m/s2) and α is surface slope derived from GPS data. We assume a constant strain rate $\dot {\varepsilon }_{zx}$ and relate $\dot {\varepsilon }_{zx}$ to the net angle of rotation ψ of the tilt sensor through the horizontal down-flow displacement ℓ
From Fig. 12, we trigonometrically calculate ψ using ℓ and sensor length L (0.205 m)
We then calculate $\dot {\varepsilon }$ from ψ