Full-physics 3-D heterogeneous simulations of electromagnetic induction fields on level and deformed sea ice

Abstract

In this paper we explore simulated responses of electromagnetic (EM) signals relative to in situ field surveys and quantify the effects that different values of conductivity in sea ice have on the EM fields. We compute EM responses of ice types with a three-dimensional (3-D) finite-volume discretization of Maxwell’s equations and present 2-D sliced visualizations of their associated EM fields at discrete frequencies. Several interesting observations result: First, since the simulator computes the fields everywhere, each gridcell acts as a receiver within the model volume, and captures the complete, coupled interactions between air, snow, sea ice and sea water as a function of their conductivity; second, visualizations demonstrate how 1-D approximations near deformed ice features are violated. But the most important new finding is that changes in conductivity affect EM field response by modifying the magnitude and spatial patterns (i.e. footprint size and shape) of current density and magnetic fields. These effects are demonstrated through a visual feature we define as ‘null lines’. Null line shape is affected by changes in conductivity near material boundaries as well as transmitter location. Our results encourage the use of null lines as a planning tool for better ground-truth field measurements near deformed ice types.

1. Introduction

Sea-ice area, thickness and volume are increasingly recognized as key variables for short-term weather and long-term climate variations in both scientific literature and public discourse (Francis, 2013). With sea ice spanning thousands of square kilometers, but only varying by meters in height, the accuracy of thickness is critical in determining sea-ice volume and monitoring its ongoing changes. To demonstrate the impact of thickness uncertainties given this large aspect ratio, an average ice thickness of 2 ± 1 m yields a range of uncertainties up to ± 5 0% in volume calculations.

Thickness measurements are collected at different scales and from different observation platforms located underwater, in and on the ice, airborne and from space, with very little integration or coincidence between platforms (Geiger and others, 2015a). The matter of thickness accuracy is an ongoing research problem intimately connected with the development of reliable techniques for rapid reconnaissance. These technical factors are being explored by posing important science questions related to changing ice patterns and their effects on the planet’s thermal stability. In this larger context, the matter of determining sea-ice thickness has become one of the priorities in Earth System observations (OSTP, 2014). One example is the need for increased accuracy of sea-ice thickness for decision-making regarding civil infrastructure of coastal polar communities (Karl and others, 2009, pp. 141-143).

The determination of sea-ice thickness using electromagnetic (EM) techniques is closely coupled to the presence of underlying conductive sea water. Several geophysical techniques are used for imaging the subsurface ice/water interface. Such methodologies include seismic dispersion within the ice (e.g. Marsan and others, 2012), electrical conductivity between the ice and the underlying water (Pfaffhuber and others, 2012), or some combination of both (Haas, 1997). Here we focus on methods based on EM induction.

To summarize, an EM induction instrument (Fig. 1) is based on magnetic dipole principles using a transmitter/ receiver pair separated by a distance r, known as the coil separation. The transmitter is coiled wire configured to generate a sinusoidal-varying electric current (Wightman and others, 2003) at a specific low frequency f (in the kHz range). The electric current sets up a time-harmonic primary magnetic field P which emits spherically from the transmitter. Nearby conductive materials respond to the emitted primary magnetic field by generating electric eddy currents which in turn create their own secondary magnetic fields S. Since the secondary magnetic fields are passively responding to the primary field, the response is called an induction or induced magnetic field (Fitterman and Labson, 2005, p. 303). For completeness, the transmitter and receiver use complex signals both containing in-phase (real) and quadrature (imaginary) components, with the secondary electric and magnetic responses being 90° out of phase with the primary fields (Wightman and others, 2003). Since the secondary field is 90° out of phase with the primary field (McNeill, 1980) in frequency-domain EM induction instruments that operate in low induction numbers (see Appendix and also the appendix of McNeill, 1980), the coil separation r is small compared to the skin depth of the surrounding materials. Essentially, the receiver coil of a low-induction number magnetic dipole pair behaves like an antenna to detect the magnetic field produced by nearby geophysical materials from eddy currents induced by the transmitter.

Fig. 1. Schematic of typical EM induction model through multiple level materials. Note that this schematic represents a vertical dipole configuration. The secondary field is induced in the receiving coil by the eddy currents that were created (induced) in each material by the transmission of a primary field from the transmitter coil, separated by a fixed length r from the receiver coil. M represents the number of layers, h is the height of the layers, and σ is the conductivity of a particular material layer.

In sea-ice applications, an EM instrument is traditionally placed either on or above the snow and/or ice surface at a height h. EM readings are traditionally calibrated to the depth of the ice/ocean interface (e.g. Geiger and others, 2015a) through the assumption that the underlying sea water is the primary conductive source, with air, snow and sea-ice conductivity assumed negligible or very small when compared to sea-water conductivity levels (Pfaffling and others, 2007).

For commercial instruments, such as the Geonics EM31-MK2 conductivity meter (hereafter referred to as EM31), the value reported to the user is known as ‘apparent conductivity’, which is proportional to the ratio of the secondary (quadrature) magnetic field divided by the primary magnetic field (McNeill, 1980). Apparent conductivity is intended to represent the conductivity a of an equivalent homogeneous earth given by the integrated contribution from all materials that are sensed by the receiver. The measurement of apparent conductivity for low induction numbers is detailed in the Appendix with parameters defined in Section 3, following McNeill (1980).

An alternative to field calibration is a layered-earth model (e.g. Wait, 1962) where the mathematical expression (see Appendix) is solved directly through the assumption of negligible conductivity for the air, snow and ice layers and one known conductivity layer (∼2 Sm^–1) for sea water. These semi-analytic formulae for layered-earth levels are simplified by an approximation through a series of exponential functions for distance between the antenna and the ice/ water interface (Pfaffling and others, 2007). Two critical assumptions for this approximation are (1) all layers are level and (2) conductivity values are uniform for the material within each layer. The thickness solution is rendered through digital filter techniques (e.g. Anderson, 1979) based on Hankel transforms (Balanis, 1989, appendix IV).

Because of the assumptions of uniform sub-layers (Fig. 1), only level ice grown thermodynamically in the absence of deformation meets the criteria. As soon as major physiographic features arise (e.g. ice rafting, ridges and cracks), the two critical assumptions are violated and the approximation previously discussed develops notable errors. The errors are often reported as a one-sided bias, with the deformed ice being systematically interpreted with less thickness of order 50% relative to validation drillholes (Reid and others, 2006; Pfaffling and others, 2007). These errors are a significant problem because true thickness is not accurately reported and thickness distribution is not well represented, though the thickness distribution does conserve volume. Geiger and others (2015b) examine this problem and demonstrate that the estimates of deformed ice are smoothed in thick locations while simultaneously measured as thicker retrievals in surrounding thinner ice as an anticipated result of averaging and an essential condition to conserve volume. Since these errors are reported as a bias (i.e. a one-sided uncertainty), they can also be classified as having two onesided biases of different magnitudes and therefore exhibit a skewed uncertainty distribution across any thickness profile. For clarity, a bias is defined as the system error between the mean value and the assigned reference value, and uncertainty is the estimate of the error for a given experiment (White, 2008). Without a point-by-point tag identifying when a measurement is biased positive or negative, both bias directions are incorporated into a common uncertainty that is much larger than that of level-ice uncertainty located a sufficient distance away from deformed features. For users of collected EM data, the different forms of error reporting (bias vs uncertainty, one-sided vs two-sided, skewed vs not), varying instrument footprint sizes, and non-standardized sampling intervals create a great deal of confusion. In particular, most of the above error parameters are not included in a corresponding data archive. Instead, the user is asked to see peer-reviewed journal articles that contain executive summaries of methods but lack the depth (and often clarity) needed to evaluate the data quality beyond a general bulk uncertainty.

Given the matters just summarized, this paper seeks to understand how physiographic features affect instrument response. As a first look into this problem, we use a three-dimensional (3-D) EM induction model that renders the steady-state heterogeneous responses through Maxwell’s equations (Weiss, 2013). Since airborne EM induction systems are already advancing with 3-D detection capabilities (Pfaffhuber and others, 2012), we focus on the simulated responses of EM signals relative to in situ field surveys to improve the ground-truth capability in support of airborne, spaceborne and underwater remote-sensing instruments.

2. Numerical Experiment

The EM response of a generalized 3-D air/snow/ice/ocean system is modeled here using a finite-volume heterogeneous solution to Maxwell’s equations (Weiss, 2013). The model is called Project APhiD, where APhiD stands for magnetic vector potential A and electric scalar ϕ (‘Phi’) Decomposition, which fully accounts for both ohmic conduction (σ) and displacement current (J) effects (Weiss, 2013). The model domain is composed of a rectilinear grid, the cells of which are each endowed with a particular conductivity (σ) / permittivity (ɛ) pair of material properties. Magnetic permeability (μ) variations are expected to be negligible in sea-ice/water systems and thereby safely ignored (Weiss, 2013).

Model cells with defined electric properties are used to describe each material within the model domain. Specific model inputs are: (1) definition of grid geometry (cell shape and size) and (2) electric properties of each gridcell. A Cartesian Yee gridcell (Yee, 1966) describes the location and direction of EM field components (Fig. 2). The transmitter is defined as a superposition of electric current elements lying along cell edges. In this way, a loop antenna, like that in the EM31 instrument, is approximated in the model by a set of four points along a chosen cell face within the model domain. With the source and conductivity properties at each gridcell thus defined, EM potentials throughout the model domain are computed by inverting the finite-volume system of equations using a matrix-free iterative scheme which is subsequently differentiated to yield estimates of the observable fields (Weiss, 2013). Because the simulator computes the fields everywhere, each gridcell behaves like a receiver within the model volume, with each receiver point being a distance r (Fig. 1) from the transmitter. From this perspective, the response of many receivers relative to one transmitter can be visualized from each steady-state solution.

Fig. 2. Yee gridcell configuration for APhiD. H and E represent the magnetic and electric field intensities (Am^–1 and Vm^–1), and their direction. The black and blue dots indicate nodes where grid corners and faces align. Indices i, j and k relate Cartesian space along x, y and z directions, respectively.

To apply APhiD to a sea-ice environment (Fig. 3), the input parameters are as follows. First, a source is needed. An inductive source at a specific frequency is chosen with properties matching the EM31, which has an operating frequency of 9.8 kHz (Eicken and others, 2001). Model simulation runs in this study are configured with the transmitter positioned 0.5 m above the surface of the snow to mimic real carry heights in the field. Second, a grid of certain dimensions needs to be constructed to determine the field response of the instrument. In this study, a 100 x 100 x 150 cell grid is set up with resolution of uniform size 0.5 m cubed-cells for 100 vertical cells attributed with snow, ice or water properties, and 50 vertical air cells. For simplicity, the air/ice boundary is level across the entire horizontal surface a t z = 0m. We run four simulation cases (Table 1).

Fig. 3. Model geometry. (a) The level-ice case configured to match the solution (Fig. 1) but expanded to a 3-D grid. (b) A deformed-ice case involving a simple triangular ridge below sea level (no surface deformation). The z direction is positive downward, opposite of that shown in Figure 1. H_air and H_earth show the direction of the magnetic field pulse (Am^{– 1} ) through the respective media. Transmitting coil shows the direction in which we position the transmitter at certain locations along the ridge during sequential model runs, with the center of the coil shown by the black dot. Receiving coil is each gridcell in the model volume, with the receiver coil located in the center of each cell and being a distance r from the transmitter.

Table 1. Simulation types and material properties

Two sets of simulations are initially run as controls to provide a reference relative to the existing literature for well-known level sea-ice situations. The first control run contains air and sea-water layers only. The second control run includes the previous two layers with 0.5 m of snow added above a 3 m layer of level ice in a configuration matching traditional level-earth models. The transmitter for these control runs is located in the center of the horizontal face of a gridcell just above the surface (at z = 0.5 m; Table 1).

To simulate a ridge in APhiD, a mask file, which describes the dimensions of the ridge as well as other material layers, is created based on two criteria. First, the size of the mask is determined by ΔX ₁ , ΔX₂ and ΔX₃ (Fig. 3b), where ∆X1 is the distance from the left edge of the model to the start of the ridge (i.e. 15 m across given 0.5 m resolution of each cell), ΔX² is the width of the ridge (20 m) and ΔX³ is the distance from the right side of the ridge to the right edge of the model (15 m). Second, an integer value is assigned to each gridcell to discriminate which material corresponds to that cell, such that a mask value of ‘0’ corresponds to a cell with air, ‘1’ is a cell with ice, ‘2’ is snow and ‘3’ is air. Each mask value is subsequently matched with appropriate conductivity values.

We represent an ice ridge by adding a simple triangular shape, without surface deformation, to the bottom of the ice layer such that the apex of the ridge is in the center of the mask file. The simulations are run for each ridge scenario to evaluate responses to different positions of the transmitter relative to the ridge. For the first position run in each ridge scenario, the transmitter is located halfway between the start of the model domain and the beginning of the ridge. For the second position run, the transmitter is located at the beginning of the ridge where it just starts to deepen below the level ice. For the final position run, the transmitter is located above the apex of the ridge. For simplicity, all simulations reported herein are in the vertical dipole configuration (Fig. 1), meaning that the magnetic dipole aligns vertically for coils in the horizontal plane, but the horizontal dipole configuration is also easily adapted.

3. Simulation Characteristics

Using the control runs listed (Table 1), we visualize four specific features (Fig. 4) using 2-D profile slices to represent the entire volume under study. These four specific features begin with the two well-known EM parameters (1) electric current density lines and (2) magnetic flux density lines. Additionally, we define two parameters to communicate visually interesting characteristics in the field. These two features are (3) color map of magnetic flux density quadrature component in the vertical direction (B_Z*), and (4) effects we call ‘null lines’ which manifest as concentrated gradients in the color map of magnetic-flux density quadrature each time the polarity of the pulse changes.

Fig. 4. Control runs showing simulated field responses from APhiD to a resolution of 0.5 m. Representative slices of the 50 x 50 x 75 m³ volume provided along the vertical x-z slice and tilted horizontal x-y slice between induction instrument transmitter source and (a) an air and water layer and (b) a layer of air, 0.5 m of snow, 3 m of flat level sea ice, and sea water. In the tilted horizontal x-y slice, the white semicircles are a half-cut representation of the normalized real component of the electric current density (J) in the media as induced by the transmitter magnetic field. The yellow curves emanating from the origin represent the normalized imaginary component of the magnetic flux density (B). The color map of the quadrature component of the magnetic flux density in logarithmic space is shown along x-z and tilted x-y slices (axes in m). Null lines (highlighted in (a) by black arrows) are defined as polarity changes in traveling direction of the transmit signal. These null lines indicate the shape of the magnetic field interaction into the (a) water, and (b) ice then water. Kinks in the null lines and magnetic flux density lines indicate a material discontinuity. In (a) and (c), L represents the extent/width (spot size) of the magnetic field at the points of polarity reversal. In (b), layers of snow, sea ice and sea water are indicated for clarity. (c) Close-up of the boxed region of interest from (b) to emphasize kinks, skin depth and footprint size.

3.1. Electric current density lines

By Faraday’s law, a time-varying magnetic field induces an electromotive force, which produces an electric current density in a medium. White semicircles in the tilted x-y horizontal slice (Fig. 4) represent the normalized real component of these current flowlines. Mathematically, the electric current density lines, denoted symbolically as J with units of A m^–2are related to the conductivity a in S m^–1 and electric field intensity E by Ohm’s law (Fig. 2)

(1)

3.2. Magnetic flux density lines

The magnetic flux density (and for our purposes, the normalized imaginary component of the magnetic flux density) is represented by yellow curves emanating from the origin in the x-z plane (Fig. 4a). Generally, magnetic flux density is expressed mathematically as

(2)

Here H is the magnetic field intensity (Fig. 2) and μ is the magnetic permeability, which we assume to be the permeability of free space (4π X 10^–7Hm^–1) . Assuming a source-free region (Balanis, 1989, pp. 1-31), the magnetic flux density is also related to the electric field intensity, and, in turn, electric current density, by way of Faraday’s law and Ohm’s law (Eqn (1)), where Faraday’s law is expressed as

(3)

where ω is the angular frequency (rads–¹) and is equal to 2πf, and j = / 1 which represents the imaginary component. In a physical sense, the electric current density lines and magnetic flux density lines are related to each other since their respective fields are transverse waves, where they are mutually perpendicular and also perpendicular to the direction of propagation (Witten, 2006, pp. 139-146 and 234-236). Additionally, magnetic flux density B (Wbm^–2 or T) has in-phase and quadrature components of

(4)

where * indicates the quadrature (imaginary) components produced by the induced eddy currents, and terms without asterisk are the in-phase or real components produced by the transmitter (McNeill, 1980; Lyons, 2008).The strength of the quadrature component of the magnetic flux density in the z-direction (i.e. B*_z) is scaled on the color map (Fig. 4a), with red indicating strongest magnetic flux density. The values of the color map are stated mathematically as

(5)

It is also important to note that the higher values of magnetic flux density occur at the location of the strongest conductivity response (see Appendix).

3.3. Color map

As stated previously, the color map (Fig. 4) represents the strength of the magnetic flux density’s quadrature component in the z-direction (C ). A change in color shows the exponential decay of the transmitter pulse as it travels through the sea ice and sea water.

3.4. Null lines

When B_z* on the right-hand side of Eqn (5) equals zero (i.e. C= log₁₀|0|), the magnetic field changes polarity during positive and negative cycles as the alternating-current transmitter signal travels in the downward (positive) z-direction. We refer to these polarity changes as ‘null lines’ because the absolute value of an oscillating signal creates strong gradients, which manifest as strong horizontal lines from a vertically transmitting oscillating pulse of the magnetic field. For illustration purposes, we provide a one-dimensional (1-D) heuristic example (Fig. 5) in the form of the well-known damped oscillator of an alternating normalized current sine wave (i.e. what the EM31 generates from the transmitter coil). Assuming a 1-D travelling wave propagating downward in the x-z plane, a typical decaying (or attenuating) travelling wave within a lossy medium (e.g. sea ice and sea water) has the general solution form of

(6)

Fig. 5. Decaying sinusoidal wave. (a) A typical sinusoidal wave as it decays when penetrating through a material. Null lines occur at polarity reversals, as denoted by the red circles on the zero line of the amplitude. (b) Absolute value of the decaying waveform in (a), and also denotes where the null lines/polarity reversals occur at the red circles. The positive envelope is shown for clarity. (c) Plot, in logarithmic space, of how the actual null lines occur when the log of the decaying waveform is taken. Here the blue waveform (x ₅) possesses ‘spikes’ when the logarithm approaches –∞, resulting in a null line. The logarithm of the envelopes (Fig. 5a) forms a straight line of the maximums of the logarithm of the decaying waveform.

where A (m) is the amplitude, ∆₀ (m) is the initial wave amplitude, α (Npm^–1) is the attenuation constant, z is the direction the wave is travelling in, t (s) is the time, β (rad m^–1) is the phase constant and ϕ (rad) is the reference phase. It should be noted that

(7)

where A (m) is the spatial wavelength of the wave. We choose sine dependence to express the imaginary component as in Euler’s identity,

(8)

where θ refers to the phase angle (Ulaby, 2004, p.22).

For our study, this general form reduces to the steady-state solution for a z-directed wave, with A0 = 0.5m, a = 5Npm–¹ f=10Hz, λ = 0.2m, t=0s and ϕ= ₀ = 0rad. Hence, the waveform and its respective envelopes are described mathematically as

(9)

(10)

(11)

where the envelopes show the decrease in amplitude with distance (Witten, 2006, p. 133).

As the wave travels through a material, the signal attenuates, but still maintains the waveform. When we take the absolute value of the waveform of Eqn (5), we see that the waveform shape is no longer represented in negative space (i.e. retaining only positive values; Fig. 5b). When the waveform intersects zero amplitude, the polarity changes. When the logarithm of zero is taken as in Eqn (5), we start to notice strong gradient lines asymptotically approaching –∞. In equation form (Fig. 5c), we state mathematically that

(12)

When we take the absolute value of either envelope in Eqn (10) or (11), and then take the logarithm, we note that the result forms a straight line of the maximums of the logarithm of the decaying waveform. This equation is expressed as

(13)

(14)

for the envelope in Eqns (10) and (11). Note that with these equations, ‘spikes’ occur each time the waveform amplitude (Fig. 5a) equals zero, which indicate a change in polarity; thus, a null line is defined. From a geophysical perspective, null lines appear near material boundaries where transmitted signals refract off the sharp material gradient interface.

Another important relationship to describe with the characteristics of this model is skin depth. Skin depth describes the effective penetration depth of an emitted signal through each material, and is expressed as

(15)

It is important to note here that skin depth is not only a vertical penetration parameter, but also a measure of horizontal penetration such that skin depth means the penetration depth through a material in any direction. The ability of APhiD to resolve 3-D structure makes it possible to explore horizontal issues that traditional 1-D level-earth models could not (Fig. 1). Most importantly, the mapping of skin-depth patterns through the location and visualization of null lines provides us with an effective tool to characterize material conductivities, which are presented next.

4. Simulation Results

The first control run is of the air/water interface only (Fig. 4a). Here we see null lines forming a kink at the air/ water interface. We subsequently label this feature a ‘kink’ because of the slight bend in the null line shape at the material discontinuity between the air and water layers. Another important shape to these null lines is their width relative to the coil spacing r. Null line width L, also referred to as footprint (term derived from footprint discussed in Geiger and others, 2015a) or spot size from beam optics (Saleh and Teich, 1991, p. 85), defines the footprint size of each polarity change at each depth. Comparing L (Fig. 4a and b), footprint is sensitive to each material’s conductivity, which is important for showing the relationship between material conductivity and the scale of the EM field response to different materials.

Additionally, for the air/sea-ice/water solution (Fig. 4b), the magnetic flux density lines also show the kink at the ice/ water interface, 3 m below sea level. Notice that the kink here is at a sharper angle in the air/ice/water case than the air/water case (Fig. 4a). Accordingly, this simulation shows that the footprint L of the magnetic field widens with the inclusion of a layer of level sea ice, compared with case for no sea ice. This result is something unanticipated in prior 1-D level-earth models.

When we add a simple triangular ridge with the same conductivity as the ice in the control run (Fig. 6), the resulting EM field pattern gets complicated. First, depending on the location of the transmitter, the current density lines (Fig. 6; white lines, normalized real component of J) outline the ridge structure. Meanwhile, null lines are skewed in various directions depending on the location of the transmitter, with kinks appearing within homogeneous layers. Furthermore, results (Fig. 7) show variations due to differences in sea-ice conductivity for two different ridge cases. The level and ridged ice have two different conductivities (Fig. 7) following the observations of Pfaffhuber and others (2012). These differences impact the current density lines such that they are flatter in shape than the previous simulation (Fig. 6), and the null lines take on a different shape. Additionally, when the transmitter is directly over the apex of the ridge, the footprint (L) is narrower than in the previous results (Fig. 6).

Fig. 6. (a–f) Simulated current density lines (normalized real component), magnetic flux density lines (normalized imaginary component) and color map of the quadrature component of the magnetic flux density in logarithmic space for multilayer structure with a MY ice ridge, using Haas and Druckenmiller (2009) values (0.020 S m^–1 for both level and ridged ice). This configuration has air, 0.5 m of snow, 3 m of level sea ice, a MY ridge, and sea water, scaled (Fig. 4) with the transmitter loop as a black ‘source’ dot shown in an approximate horizontal location for clarity. In (c), the layers of snow, sea ice, sea water and an ice ridge are labeled for clarity. Field line results are shown from two perspectives with source at three locations. Properties of simulation described in Table 1, listed as simulation No. 3.

Fig. 7. Unconsolidated first-year deformed ridge using conductivity values (0.170 S m^{– 1} for flat ice, and 0.5 S m^{– 1} for ridged ice) from Pfaffhuber and others (2012), following Figures 4 and 6 (simulation No. 4 in Table 1). (a–c) Front view of the EM field mapping of a ridge structure with different conductivity. (d–f) Perspective view from below the sea-water/ice interface. Note that the magnetic flux density lines (normalized imaginary component) are more compressed relative to those in Figure 6 since the ice is more conductive in this scenario.

As a direct comparison of the various simulation outcomes, we plot the ridge effect (RE) differences of the imaginary component of vertical B (Fig. 8), i.e. B*z, at each gridcell between the level-ice control case (Fig. 4b) and both ridged cases (Fig. 6 and 7), and cast the absolute value of the difference into logarithmic form. In a mathematical sense,

(16)

Fig. 8. Ridge effect differences in null line shape due to the difference between the imaginary components of vertical B (B _z *) cast into logarithmic form. Results of imaginary vertical B component (B _z *) from ridged ice cases (Fig. 6 and 7) are subtracted from the imaginary vertical B component (B*) of the control run (Fig. 4b), then cast the absolute value of the difference into logarithmic form to demonstrate change in structure of the EM field lines between simulations. (a-c) The results from Figure 6 subtracted from those of Figure 4b. (d-f) The results from Figure 7 subtracted from those in Figure 4b. Transmitter is shown with red outlined dot and placed in its approximate horizontal location for clarity. Note that the strongest differences occur in (b) and (e) with the transmitter at the edge of the ridge rather than the apex.

The impact of a ridge is more apparent in the consolidated multi-year (MY) ice ridge case (Fig. 8a–c) since the ice is less conductive than in the unconsolidated first-year ridge case. Interestingly, the null lines change as the position of the transmitter moves, relative to the deformed-ice position, across the x-axis.

5. Discussion

Low-frequency induction results examined here identify the ice/water interface through the exponential decay of long-wavelength secondary-eddy-field responses in the near field, i.e. distances much less than a wavelength. This exponential decay generates a system of eddy currents induced within each material layer. The result is a sensitive measure of distance between the receiving antenna and material conductivity along any direction. As such, the sensitivity of eddy-field responses is comparable to ground-penetrating radar (GPR) responses since both concepts are based on the detection of a signal through a medium from a transmitting device (Witten, 2006, pp. 166–170 and 234–236). However, when compared to EM induction instruments, GPR has difficulties with high-salinity ice (brine-inclusive ice) as it will attenuate the signal and cause scattering (Haas and others, 2003, pp. 18–19; Pfaffling, 2006, p. 15). In line with this thought, it seems that EM induction, through the results of these simulations, may also be affected from scattering, perhaps more so than first thought. Hence, this study is only beginning to simulate a number of new and interesting responses of EM systems when explored as 3-D responses.

Some of the positive outcomes of the model results are as follows. Firstly, the modeling results shown herein can be used to plan more effective field experiments before making expensive excursions to the Arctic; specifically seasonal and regional sensitivities related to strong conductivity gradients as well as site selection, line survey selection methods, and approaches to physiographic features. In particular, model simulators such as APhiD provide a capability to numerically test novel sensor packages that measure the complete, three-component induction field for a given transmitter antenna. These novel packages can be designed either through a range of fixed frequencies, or preferably as a transient pulse. As precedent, we note that such an approach was adopted in the early 2000s by service providers in the hydrocarbon well-logging industry, such as Schlumberger, Baker Atlas and Halliburton (Anderson and Barber, 1995; Beard and others, 1996; Beste and others, 2000), when it became apparent that geological anisotropy confounded traditional methods of well-log data analysis. A decade later, deployment of three-component, multi-sensor induction logging tools is the state of practice in difficult drilling environments, geo-steering and measurement-while-drilling; albeit still proprietary. We therefore see the results of this study as potentially leading to a similar shift in sea-ice measurement technologies, thus inviting future mapping work and analysis on the full 3-D nature of climate/sea-ice dynamics.

Secondly, airborne 3-D modeling work in sea ice was attempted before (e.g. Liu and Becker, 1990), but, after extensive and advanced computing, only minor improvements were made in analyzingfield data (Pfaffling and others, 2007). Since computers have increased in processing power overtime, a new capability now exists with APhiD to produce 3-D EM models. Like the results here, the output produces interesting results that not only show the current density lines and the magnetic flux density lines, but also show how the entire field reacts as a function of 3-D distributed material conductivities. As shown in the Appendix, the relationship between the model output and material conductivities provides an opportunity to increase understanding between geophysical properties (e.g. material conductivity) and instrument responses (e.g. apparent conductivity). When the field encounters level sea ice instead of an air/water interface (Fig. 4), the footprint size increases along the horizontal. This footprint is important because the instruments used in the field assume a footprint based primarily on the carrying height of the instrument (Geiger and others, 2015a). Conversely, in these findings, we see that the footprint of a pulse varies considerably based on the conductivity of the material, so the footprint of the instrument is going to be sensitive to ice type, season of year, temperature and other environmental variables not currently formulated in any EM models. Further studies are underway to quantify how L varies as a function of sea-ice conductivity, with the intention to develop an algorithm characterizing these changes from parameterized physical properties.

Thirdly, there is the potential use of APhiD as a planning tool for positioning instruments in the vicinity of sea-ice ridges to improve ground-truth data collection best practices. There are considerable patterns (Fig. 6 and 7) that are reminiscent of refraction and physiographic interference at the beginning of the ridge. We can postulate that the physiographic feature, i.e. the ice ridge, induces both refraction and interference patterns depending on the ridge shape and conductivity of the ice. To explain, both of these patterns are based on the Huygens-Fresnel principle. Refraction is where the EM waves ‘bend’ – visually similar to the appearance of a pencil in a glass of water - due to the wave speed change across the boundary of two different media dependent on its material properties, such as conductivity (Witten, 2006, pp. 215-216). Interference is the result of the spherical wave generated by the wavefront – in our case, the wavefront is the transmitter coil – where the envelope of the spherical waves constitutes a new wavefront via superposition (Saleh and Teich, 1991, p. 121). The only time that these patterns do not occur is when the transmitter is above the ridge apex, in which case ridge symmetry matches the field shape waveform symmetry.

Hence, the beginning of the ridge may be a promising area for an EM refractive and interference process study, while the peak of the ridge is more important to EM calibration. APhiD also provides a means to rethink instrument designs for in situ measurements to capture 3-D data by leveraging understanding about changing footprints (∆L) as a function of sea-ice conductivities (or lack of conductivity in the associated snow layer, which also impacts footprint size). The magnetic flux density lines bend away from the vertical when a pressure ridge is present since the ridge is less conductive then the surrounding sea water. The MY ice (Fig. 6) is less conductive than the unconsolidated ridges (Fig. 7), hence the magnetic flux density lines are further away from the vertical in Figure 6. The tendency for these lines to appear further away provides potential insight when interpreting airborne EM responses due to spatial resolution such as NASA’s IceBridge program (Gardner and others, 2012) and the Multi-sensor Airborne Sea Ice Explorer (MAiSIE; Pfaffhuber and others, 2012). A benefit of using APhiD is that it is a tool that can be modified for any number of scenarios. Further work with APhiD includes the integrations of material conductivity for direct comparison with apparent conductivity both modeled and in the field.

6. Conclusions

A set of simulations of the EM response of various sea-ice types has been discussed here in order to determine how an EM induction instrument, commonly used in sea-ice field studies, interacts with different conductivity values within different media. Current approximations using the 1-D layered-earth solution assume that material layers above the conductive sea water are essentially negligible in terms of conductivity. However, with the results presented in this paper, we can conclude that changes in conductivity, especially at the sea-ice/sea-water interface, do have an effect on the overall EM field response. A simple change of the location of the transmitter, with respect to deformed ice, also leads to interesting EM field responses as demonstrated in the visualizations of the deformed-ice simulation cases.

In summation, we identified four key findings. First, using a full-physics, heterogeneous finite-volume EM model, we demonstrate how the inherent assumptions in existing 1-D model approximations are violated when physiographic features are present. Second, APhiD-modeled EM parameters can be combined to show where null lines are located. Third, kinks in the null lines not only develop at material interfaces, they also develop uniquely through field patterns when physiographic features are present. While these patterns are more pronounced when they encounter the beginning of a ridge as in the simulation cases described above, further study is warranted to explore if these patterns can aid in determining the shape of the ridge in addition to the thickness. Finally, the most important outcome is that sea-ice conductivity has an important role to play in the horizontal extent of EM footprint sizes and therefore is a key parameter for interpreting EM thickness retrievals from field campaigns. These simulation results can be explored further to develop new in situ instruments, improve ground-truth calibration and validate airborne, spaceborne and underwater instruments.