A Network Model for Electrical Transport in Sea Ice

J. Zhu, K. M. Golden, A. Gully, and C. Sampson

University of Utah, Department of Mathematics 155 S 1400 E RM 233, Salt Lake City, UT, 84112-0090 USA

Abstract

We develop a network model for the electrical conductivity of sea ice, and obtain close agreement with the results of experiments we conducted in Antarctica. Monitoring the thickness of sea ice is an important tool in assessing the impact of global warming on Earth’s polar regions, and most methods of measuring ice thickness depend on detailed knowledge of its electrical properties.

Key words: sea ice, electrical conductivity, network model, sea ice thickness PACS Numbers: 92.10.Rw, 77.84.Lf, 77.22.Ch, 72.80.Tm, 93.85.Jk, 93.30.Li, 93.30.Ca

1. Introduction

Sea ice is a sensitive indicator of climate change, as well as a critical component of Earth’s climate system. Determining the thickness distribution of the polar sea ice packs is a central problem in monitoring the impact of global warming. However, there is significant uncertainty in our knowledge of the ice thickness distribution and how it is changing. Not only does this uncertainty affect assessments of how the changing climate is impacting the polar regions, but it also affects predictions of global climate models, where accurate knowledge of sea ice initial conditions is essential for long term simulations.

Most methods of measuring sea ice thickness, and interpretation of the data obtained, depend on detailed knowledge of the electrical properties of the ice. Since sea ice is a composite of pure ice with brine inclusions [20,3], whose volume fraction and geometry depend strongly on temperature, understanding its electrical properties is a challenging problem in the theory of inhomogeneous materials. Here we develop a network model for the electrical conductivity of sea ice, and compare the results with direct measurements of the vertical conductivity of first year

Email addresses: zhu@math.utah.edu (J. Zhu),

golden@math.utah.edu (K. M. Golden),
gully@math.utah.edu (A. Gully),
christian.sampson@gmail.com (C. Sampson).

sea ice we made during the 2007 Sea Ice Physics and Ecosystem eXperiment (SIPEX) expedition off the coast of East Antarctica, from the Australian icebreaker Aurora Australis.

Early DC resistivity measurements on sea ice were aimed at determining ice thickness [5,18,19]. Initially all these studies employed surface soundings using 4 electrodes in either the Wenner or Schlumberger configurations, although Timco [19] later used vertically arranged electrodes in the side of an ice pit. Later measurements in the Antarctic were also reported [2]. The anisotropic nature of the resistivity of sea ice, due to the preferential vertical alignment and connectivity of brine pores, leads to such measurements significantly underestimating the ice thickness.

More promising determinations of sea ice thickness have been achieved using low frequency electromagnetic (EM) techniques [14,11,13,21,17]. The technique relies on a time varying primary magnetic field (generated by a transmitter coil) inducing eddy currents in the seawater beneath the comparatively resistive ice. The secondary magnetic field produced is sensed by a receiver coil, determining an apparent conductivity which results essentially from an integration over the vertical distance between the instrument and induced currents. The thickness is found using empirical relationships [12], with good results for smooth ice and underestimates near ridges [12]. The technique is adaptable to continuous measurements being made either from a helicopter or ship [11].

Previous measurements of the conductivity of sea ice relied almost exclusively on indirect methods which mix the horizontal and vertical components. Moreover, these indirect means make it difficult to accurately recover the dependence of the conductivity on the properties of the brine microstructure, namely, its brine volume fraction ϕ, which depends on the temperature T and salinity S of the ice [4,20,3]. During the 2007 Australian SIPEX expedition, Golden and Gully extracted cylindrical cores of sea ice and made vertical conductivity measurements along these cores using metal probes attached to a Yew Earth Resistance Tester, as described in [9]. We also measured salinity and temperature along each core in order to relate the electrical measurements to microstructural data [15,7,8,16].

Part of our motivation for focusing on the vertical component of the electrical conductivity is that it is closely related to the vertical component of the fluid permeability of sea ice. Fluid transport in sea ice mediates a broad range of processes such as the growth and decay of seasonal ice, the evolution of melt ponds which determine ice pack albedo, and bio-mass build-up [8,6]. Our work here will help lay the foundation for electrically monitoring fluid transport in sea ice. In fact, the random resistor network model we develop here is based on the random pipe network we used previously to model the fluid permeability of sea ice [22]. Statistical information about the brine microstructure [15,7,8,16] is used to determine the statistical distributions of the resistors in the electrical network.

2. The Network Model for the Effective Conductivity of Sea Ice

In this model, we consider a piece of sea ice with an averaged brine volume fraction ϕ, and focus on the effect of brine structure on the electric conduction of the material. The conducting property of a medium can be summarized by its conductivity. Let Φ be the electric potential, and σ the local conductivity tensor, which may depend on the local volume fraction, temperature, and salinity. Since the current density J is related to the electric potential through J = σ∇Φ, using Kirkoff’s law and assuming the material is free of electrical charge, the equation for electrical conduction is

∇ · σ ∇Φ = 0, (1)

which is similar to the fluid pressure equation resulting from Darcy’s law

∇ · k ∇p = 0, (2)

where p is the the flow pressure and k is the permeability tensor.

Here we intend to define an effective conductivity that describes the behavior of conduction of the sea ice structure in the direction that is of importance to us. In particular, we are interested in the effective vertical conductivity defined through

J^{∗} = σ_{v} \frac{ΔΦ}{Δz} (3)

for the current density in z direction J^{∗}, and the potential difference ΔΦ over a thickness Δz.

To simulate the electric field through the conducting microstructure of sea ice, consider an ice sheet of depth D. We will model the medium in the way that is analogous to our previous work [22]. Take a thin vertical slice of horizontal thickness h and length span L. We model this ice sheet by a two dimensional lattice of nodes connected by conducting tubes. The slice has a rectangular L × D vertical cross section, which is divided into a grid with m equally spaced sections in the horizontal direction and n equally spaced sections in the vertical direction, so that L/m = D/n = h, for some large integers m and n. The model parameter h can be viewed as the dimension of a cell in which a typical brine inclusion is contained. In this network model, h will be chosen according to the sea ice we simulate, its brine volume fraction, and our computing capacities. The tubes are assumed to have circular shapes with different radii, and the current through the medium is induced by an electric potential drop Φ_{b} −Φ_{t}, where Φ_{b} and Φ_{t} are the potentials at the bottom (liquid) and the top (air) of the sea ice, with the assumption that Φ_{b} > Φ_{t} so there is an upward current flow in the medium. The cross sectional areas of the tubes chosen below generate fluid pores comparable to the brine inclusions found in young sea ice. The lattice nodes are the vertices (i, j), 0 ≤ i ≤ m, 0 ≤ j ≤ n, of a rectangular grid, as in Figure ?? (a). Nearest neighbors are connected by vertical and horizontal tubes, with a potential Φ_{i,j} defined at each node (i, j). To each node (i, j) with 0 ≤ i ≤ m − 1, 0 ≤ j ≤ n − 1, the horizontal tube to the right of (i, j) has radius R = R_{i,j}^{h}, and the vertical tube on top of (i, j) has radius R = R_{i,j}^{v}, as shown in Figure 1. Along the right edge with i = m the nodes have one vertical tube, and along the top edge the nodes have one horizontal tube (except the last).

Since the brine conductivity is substantially higher than the conductivity of the surrounding ice (on an order of about 10^{8}), we can assume that electrical conduction takes place mostly through the

Fig. 1. (a) (b)

brine tubes. The effect of conduction through pure ice will be modeled by adding a simple conducting component to the system. Unlike the permeability model, where the fluid flux depends only on the brine geometry, the electric conduction in the microstructure would include a temperature dependent local conductivity. For each tube of radius R connecting two nodes, we assume an uniform conductivity σ_{tube} for the brine, so an electric current within can be established based on the voltage drop and the cross sectional area. In our particular case we have

I = σ_{tube} A E = −σ_{tube} πR^{2} ∇Φ, (4)

where Φ is the electric potential. For each tube connecting two neighboring nodes, the potential gradient can be well approximated by the difference between the potentials at these nodes, divided by the spacing h. Given the potentials at neighboring nodes, different fluxes converging to the node (i, j) can be easily computed, and they must balance due to Kirkoff’s law. Let σ_{i,j}^{h} and σ_{i,j}^{v} denote the brine conductivity for the tubes to the right and on the top of node (i, j), respectively. This leads to the following equations,

σ_{i,j}^{v}(R_{i,j}^{v})^{2}(Φ_{i,j+1} − Φ_{i,j}) + σ_{i,j−1}^{v}(R_{i,j−1}^{v})^{2}(Φ_{i,j−1} − Φ_{i,j}) +
σ_{i,j}^{h}(R_{i,j}^{h})^{2}(Φ_{i+1,j} − Φ_{i,j}) + σ_{i−1,j}^{h}(R_{i−1,j}^{h})^{2}(Φ_{i−1,j} − Φ_{i,j}) = 0, (5)

for i = 1, …, m − 1, and j = 1, …, n − 1, with appropriate modifications on the edges of the lattice. Notice that this equation is similar to the equation derived for the fluid permeability model [22]:

(R_{i,j}^{v})^{4}(p_{i,j+1} − p_{i,j}) + (R_{i,j−1}^{v})^{4}(p_{i,j−1} − p_{i,j}) +
(R_{i,j}^{h})^{4}(p_{i+1,j} − p_{i,j}) + (R_{i−1,j}^{h})^{4}(p_{i−1,j} − p_{i,j}) = 0, (6)

here p_{i,j} is the pressure at node (i, j). We comment that the current system has coefficients dependence on the radius not as strong as that in the permeability case, but they also have an extra dependence on the local brine conductivity, which varies from one location to another due to temperature and salinity variations.

At the top of the region (j = n), we assume the boundary condition

Φ_{i,n} = Φ_{t}, (7)

and at the bottom (j = 0) of the region

Φ_{i,0} = Φ_{b}. (8)

Let I_{i,j} be the current through the vertical duct on top of the (i, j) node. The total current through the brine network system is therefore

I_{brine} = ∑_{i=0}^{m} I_{i,n−1} = π ∑_{i=0}^{m} σ_{i,n−1}^{v}(R_{i,n−1}^{v})^{2} \frac{Φ_{i,n−1} − Φ_{t}}{h} ,

(9)

where R_{i,n−1}^{v}, i = 0, …, m are the radii of the tubes connected to the top of the sea ice, and Φ_{i,n−1} are the potentials at the nodes just below the top surface. The effect of conduction through pure ice can be modeled as an additional current flow through another medium, in parallel to the brine network. The current through such a medium is

I_{pure ice} = σ_{ice} Lh(1 − β(ϕ))\frac{Φ_{t} − Φ_{b}}{D} . (10)

Here we introduce a function β(ϕ) that describes the loss of ice surface for conduction due to the brine inclusions.

When the piece of sea ice is viewed as a composite, the effective vertical conductivity σ_{v} can be defined through

J^{∗} = −σ_{v} \frac{Φ_{t} − Φ_{b}}{D} (11)

Where J^{∗} is the average current density through the medium. We make the connection between J^{∗} and the total current through

J^{∗} = \frac{I_{brine} + I_{pure ice}}{Lh} (12)

where Lh is the horizontal cross sectional area of the slice. Based on these assumptions, we have the effective conductivity

σ_{v} = \frac{πD}{Lh^{2}} ∑_{i=0}^{m} σ_{i,n−1}^{v}(R_{i,n−1}^{v})^{2} \frac{Φ_{t} − Φ_{i,n−1}}{Φ_{t} − Φ_{b}} +(1−β)σ_{ice}.

(13)

The cross section reduction form for the ice should depend on the microstructure of the sea ice under study. One possible model is to assume that half of the tubes are vertical and contributing to the vertical conduction. This leads to the following reduction factor:

β = \frac{1}{2}ϕ. (14)

The multigrid algorithm to solve the system Eq.(6) can be easily modified to solve the system Eq.(5), and the convergence is faster due to the coefficient dependence change from R^{4} to R^{2}.

3. Sea Ice Microstructure and Numerical Results

In this model, the microstructure of the sea ice slice is described by an averaged brine conductivity σ_{b}, and a collection of tubes with cross sectional area sampled from a log-normal distribution, with parameters based on measurements of brine inclusions in first year sea ice [15,?,1]. Specifically, we sample the radius R so that log A = log(πR^{2}) is normally distributed with mean μ and variance σ^{2}. We also assume that all the random radii are independent from each other. Given a particular sample of the tube radii, the volume fraction ϕ of the slice can be readily computed by

ϕ = \frac{π}{LD} ( ∑_{i=0,j=0}^{m−1,n}(R_{i,j}^{h})^{2} + ∑_{i=0,j=0}^{m,n−1}(R_{i,j}^{v})^{2}) . (15)

The goal of this study is to investigate the dependence of the effective vertical conductivity σ_{v}, and the form factor σ_{v}/σ_{b}, on the porosity ϕ, which is connected to the microstructure through Eq.(15). For consistency, it is necessary to choose the parameters μ and σ such that the desired volume fraction is arrived, and that the statistical properties of the actual sea ice are reasonably matched. To this end, we first notice that the expected value of the cross sectional area

E[A] = e^{μ+ \frac{1}{2} σ^{2}}, (16)

given our assumption about the distribution of log A. This expected value for the model should be matched to an interpolation of measured averages for the cross sectional area A as a function of brine volume fraction ϕ [10]:

⟨A⟩ = θ(ϕ) = π(7×10^{−5}+1.6×10^{−4}ϕ)^{2} m^{2}. (17)

This function approximates the dependence of the mean cross sectional area on the brine volume fraction ϕ observed by Perovich and Gow (1996) [15] in horizontal thin sections of young, primarily columnar sea ice. Since we have two parameters for the model to be determined, the matching condition leaves us with one free parameter σ, which we will choose several different values to compare with the measured sea ice data. It should be pointed out that the tomography information on the microstructure is rather insufficient and the distribution of the brine inclusions as reflected through μ and σ still requires substantial modeling. Therefore we will leave the parameter σ open for further modeling. As observed in [15], values of σ between 1 and 2 seem to be reasonable for the available topography information.

Also as observed in [?], brine channels in the sea ice become connected to a substantial level only when the brine volume fraction reaches above 5%. To reflect this behavior, we allow some randomly selected tubes to be disconnected from the system in an effort to simulate the disconnection of brine inclusions along certain directions. The capability of the multigrid algorithm makes it possible to study the general percolation problem from a new perspective. For the numerical algorithm, we choose to develope a multigrid algorithm especially targeted at this unusual linear system. Multigrid methods in general are very powerful in dealing with large linear systems resulting from elliptic equations, and they are particularly appropriate here to address different length scales involved in the microstructure. Our numerical examples demonstrate that this algorithm is robust, accurate, and efficient.

Similar to the work of fluid permeability, we will focus on the form of dependence of the effective ver-

Fig. 2. (a) (b)

Table 1 Tube Disconnection Probability

ϕ  0.025 0.05 0.075 0.1 0.0125
P rob(disconnection) 0.7 0.5 0.3 0.15  0

merical algorithm will encounter difficulties in some of the random configurations.

In Figure 2, we show the calculated values of the effective vertical conductivity according to Eq.(13), with two different choices for the brine inclusion distribution parameter σ, as compared to the measured vertical conductivity of the sea ice. Obviously the numerical results from the choice σ = 1 agree much better with the measurements than that of σ = 0.5. This seems to agree with the findings about the observed statistics of brine inclusions of Perovich and Gow [15].

A more fundamental quantity for our study is the form factor, the relative conductivity defined as σ_{b}/σ_{b}. In Figure 3, we plot the values of the form factor for our computed vertical conductivity and our measured sea ice conductivity in a log-log graph. Of major interest from the graph is an estimate of the exponent m and to see how our model agrees with the measured quantities. From this graph, with the choice of disconnect probability, we have achieved an excellent agreement for m around 2. A fine tuning of the various parameters to reflect the microstructure if more detailed tomography information is available should lead to more accurate predictions of the exponent m.

4. Conclusions

Acknowledgements: We are grateful for the support provided by the Division of Atmospheric Sciences and the Division of Mathematical Sciences at the National Science Foundation, through the Collaboration in Mathematical Geosciences (CMG) Program. Adam Gully gratefully acknowledges the support of the NSF Research Experiences for Undergraduates (REU) Program and an NSF VIGRE graduate fellowship.

References

[1] C. Bock and H. Eicken. A magnetic resonance study of temperature-dependent microstructural evolution and self-diffusion of water in Arctic first-year sea ice. Ann. Glaciol., 40:179–184, 2005.

[2] R. G. Buckley, M. P. Staines, and W. H. Robinson. In situ measurements of the resistivity of Antarctic sea ice. Cold Reg. Sci. Tech., 12:285–290, 1986.

[3] H. Eicken. Growth, microstructure and properties of sea ice. In D. N. Thomas and G. S. Dieckmann, editors, Sea Ice: An Introduction to its Physics, Chemistry, Biology and Geology, pages 22–81. Blackwell, Oxford, 2003.

[4] G. Frankenstein and R. Garner. Equations for determing the brine volume of sea ice from −0.5^{o} to −22.9^{o} C. J. Glaciol., 6(48):943–944, 1967.

[5] K. Fujino and Y. Suzuki. An attempt to estimate the thickness of sea ice by electrical resistivity method ii. Low Temp. Sci., A21:151–157, 1963.

[6] K. M. Golden. Climate change and the mathematics of transport in sea ice. Notices of the American Mathematical Society, 56(5):562–584 and issue cover, 2009.

[7] K. M. Golden, S. F. Ackley, and V. I. Lytle. The percolation phase transition in sea ice. Science, 282:2238–2241, 1998.

[8] K. M. Golden, H. Eicken, A. L. Heaton, J. Miner, D. Pringle, and J. Zhu. Thermal evolution of permeability and microstructure in sea ice. Geophys. Res. Lett., 34:L16501 (6 pages and issue cover), doi:10.1029/2007GL030447, 2007.

[9] K. M. Golden, A. Gully, C. Sampson, J. Zhu, A. P. Worby, and J. Reid. Theory and measurements of

Preprint submitted to Elsevier 28 September 2009 electrical conductivity in first year Antarctic sea ice. in preparation for a special issue of Deep Sea Res., 2009.

[10] K. M. Golden, A. L. Heaton, H. Eicken, and V. I. Lytle. Void bounds for fluid transport in sea ice. Accepted for publication in a special issue of Mechanics of Materials on Advances in Disordered Materials, 2005.

[11] C. Haas. Evaluation of ship-based electromagneticinductive thickness measurements of summer sea-ice in the Bellingshausen and Amundsen Seas, Antarctica. Cold Reg. Sci. Tech., 27:1–16, 1998.

[12] C. Haas. Dynamics versus thermodynamics: The sea ice thickness distribution. In D. N. Thomas and G. S. Dieckmann, editors, Sea Ice: An Introduction to its Physics, Chemistry, Biology and Geology, pages 82–111. Blackwell, Oxford, 2003.

[13] C. Haas. Late-summer sea ice thickness variability in the Arctic Transpolar Drift 1991-2001 derived from groundbased electromagnetic sounding. Geophys. Res. Lett., 31:L09402, doi:10.1029/2007GL030447, 2004.

[14] C. Haas, S. Gerland, H. Eicken, and H. Miller. Comparison of sea-ice thickness measurements under summer and winter conditions in the Arctic using a small electromagnetic induction device. Geophysics, 62:749757, 1997.

[15] D. K. Perovich and A. J. Gow. A quantitative description of sea ice inclusions. J. Geophys. Res., 101(C8):18,327– 18,343, 1996.

[16] D. J. Pringle, J. E. Miner, H. Eicken, and K. M. Golden. Pore-space percolation in sea ice single crystals. Journal of Geophysical Research C, in press.

[17] J. E. Reid, A. Pfaffling, A. P. Worby, and J. R. Bishop. In situ measurements of the direct-current conductivity of Antarctic sea ice: implications for airborne electromagnetic sounding of sea-ice thickness. Ann. Glaciol., 44:217–223, 2006.

[18] F. Thyssen, H. Kohnen, M. V. Cowan, and G. W. Timco. DC resistivity measurements on the sea ice near pond inlet. Polarforschung, 44:117–126, 1974.

[19] G. W. Timco. An analysis of the in-situ resistivity of sea ice in terms of its microstructure. J. Glaciol., 22:461– 471, 1979.

[20] W. F. Weeks and S. F. Ackley. The growth, structure and properties of sea ice. In N. Untersteiner, editor, The Geophysics of Sea Ice, pages 9–164. Plenum Press, New York, 1986.

[21] A. P. Worby, P. W. Griffin, V. I. Lytle, and R. A. Massom. On the use of electromagnetic induction sounding to determine winter and spring sea ice thickness in the Antarctic. Cold Reg. Sci. Tech., 29:49–58, 1999.

[22] J. Zhu, A. Jabini, K. M. Golden, H. Eicken, and M. Morris. A network model for fluid transport in sea ice. Ann. Glac., 44:129–133, 2006.