Improved 1D inversions for sea ice thickness and conductivity from electromagnetic induction data: Inclusion of nonlinearities caused by passive bucking

Along with the shrinking sea ice extent (Stroeve et al., 2012; Comiso and Hall, 2014), Arctic sea ice is thinning (Lindsay and Zhang, 2005; Haas et al., 2008; Kwok et al., 2009) and becoming younger (Maslanik et al., 2011). The change toward a younger sea ice cover leads to changes in the average internal structure. For example, first-year sea ice shows higher absorption rates of shortwave radiation than does older sea ice, leading to a stronger internal warming and melting (Perovich et al., 2011; Nicolaus et al., 2012). This in turn results in an increase of the connectivity of pores and brine channels, which governs the exchange of gases (Gosink et al., 1976) and nutrients (Krembs et al., 2011) and the mechanical strength (Kovacs, 1996). Here, we develop and present a methodology to improve the use of electromagnetic (EM) induction sounding instruments to retrieve information related to the internal structure of different sea ice types, based on the measurements of bulk electrical sea ice conductivity, hereinafter referred to as conductivity.

The conductivity and porosity of a sea ice layer are generally linked (Archie, 1942). Pores of young sea ice are usually filled with highly saline brine left from the formation process, whereas in older, multiyear sea ice, this brine is usually replaced by freshwater due to meltwater flushing (Eicken et al., 2002; Vancoppenolle et al., 2007) or gravity drainage (Niedrauer and Martin, 1979; Notz and Worster, 2008). This also means that older, multiyear sea ice is electrically more resistive, whereas young, porous sea ice still has a high conductivity due to brine inclusions. Summer sea ice often shows macroscale pores because brine channels are widened and connected by melt processes. These pores then contain a mix of meltwater and brine, which change the geochemistry of sea ice and the availability of nutrients for biological processes (Thomas and Dieckmann, 2009; Vancoppenolle et al., 2010). Mapping sea ice layer porosity by proxy measurements of sea ice conductivity would therefore enable improved process studies of sea ice mass balance and especially complex biogeochemical cycles through the atmosphere/sea ice/ocean interfaces. Information on the spatial variability, large-scale gradients, and the annual cycle of internal sea ice properties is only sparsely available (e.g., Golden et al., 2007; Pringle et al., 2009). Current methods to obtain information on sea ice porosity are based on manual coring, whereas the high spatial variability in the sea ice cover calls for a rapid and nondestructive method.

The conductivity of sea ice, and in particular its large contrast with that of the ocean, is routinely used for sea ice thickness retrieval by frequency-domain EM induction sounders, which usually use a single frequency in a range from 4 to 10 kHz. This method is deployed from helicopters (Kovacs and Holladay, 1990; Haas et al., 2009), aircraft (Haas et al., 2010), ship-mounted booms (Haas, 1998; Reid et al., 2003), or directly on the sea ice surface for high-resolution case studies (Eicken et al., 2001; Haas, 2004; Druckenmiller et al., 2009; Weissling et al., 2011). The retrieval algorithm assumes that sea ice and snow can be described as a single, flat 1D layer, which is usually the case for level sea ice, and that the effect of sea ice and snow conductivity on the EM response is negligible compared with the contribution from the underlying seawater. The EM response can be approximated with an exponential relation to the distance of the sea ice/seawater interface of either the in-phase or quadrature (out-of-phase) component, parameterized from calibration measurements or 1D forward models. The inverse of the exponential function then directly relates the EM response to the total thickness (snow and sea ice) for the given height of the sensor above the surface (Kovacs and Morey, 1991). We refer to this method as empirical exponential (EMPEX) according to Pfaffling et al. (2007). This approach has proven to be a robust methodology for empirical sea ice thickness retrieval, but it inherently has to neglect any changes of the sea ice conductivity.

The assumption of a perfectly resistive, single layer of sea ice and snow is justified in most cases because the influence of the sea ice conductivity of dry, cold, level sea ice with typical thicknesses of 1–2 m is small in the traditionally used frequency range (Pfaffling and Reid, 2009). However, there are cases in which the sea ice conductivity is desired for sea ice type classification and even becomes a dominant factor, e.g., for summer melt conditions or the so-called subice platelet layer. The latter is a layer of loose ice crystals (ice platelets) with a bulk conductivity between those of sea ice and ocean water (Rack et al., 2013; Hunkeler et al., 2015), and it may be located underneath coastal sea ice in Antarctica. Ice platelets originate from supercooled water in cavities below ice shelves and are a main contributor to the local sea ice mass balance. A mapping of this special sea ice type is still a challenge due to the scarcity of available data (Hoppmann et al., 2015; Langhorne et al., 2015). In addition, the sea ice thickness can vary on instrument subfootprint scales, and not all sea ice structures can be sufficiently described by a 1D geometry approximation within the footprint. The traditional 1D assumption leads to thickness biases for deformed sea ice, partly due to the geometric effect of subfootprint-scale thickness variability (Kovacs et al., 1995; Reid et al., 2003; Pfaffhuber et al., 2012) and partly due to the intrusion of saline seawater into the sea ice layer (Reid et al., 2006).

Operational needs for in situ high-resolution sea ice surveys require a small sensor package that can be towed by a person on foot or by snowmobile. It is a particular characteristic of such surveys that measurements are taken in a highly conductive environment given that the sensor is usually a few meters above seawater. The conductivities of sea ice layers we aim to resolve (⁠$0.01−0.2S/m$⁠) can be 1–2 orders of magnitude smaller than the conductivity of the ocean (e.g., $2.7S/m$⁠). Distances between the sensor and the saline ocean water can easily become smaller than the dimensions of the sensor itself. Previous EM sea ice studies have accounted for sea ice conductivities, such as Pfaffling and Reid (2009), who reveal that the lower frequency (3.68 kHz) of an airborne device is mostly sensitive to sea ice thickness, whereas higher frequencies (112 kHz) are needed to resolve sea ice conductivity. However, the difference in conductivity between first-year sea ice (⁠$0.05S/m$⁠) and multiyear sea ice (⁠$<0.01S/m$⁠) could not be resolved for sea ice thinner than 2 m. In contrast, Holladay et al. (1998) use airborne EM induction sounding to distinguish between first- and multiyear sea ice using 30 and 90 kHz. Reid et al. (2003) conclude that it is necessary for sea ice thickness retrieval to account for conductive sea ice.

Empirical analytic approaches for processing single-frequency EM data are insufficient to resolve thickness and conductivity of even a single layer. Therefore, we assess the prospect of using 1D inversion algorithms on multifrequency EM data, in which (1) conductivity within fixed layers and (2) conductivity and thickness of sea ice layers are adapted to minimize the difference between the signal response and forward models. Previously, Hunkeler et al. (2015) perform case studies with a ground-based device (GEM-2, Geophex Ltd.) at frequencies of 63.03 and 93.09 kHz and a coil spacing of 1.66 m, being able to resolve sea ice conductivities in agreement with estimates from sea ice cores. The assessment of sea ice conductivity in this study is based on the same small-coil instrument in which small coil refers to the transmitter (Tx) and receiver (Rx) coils that have small diameters compared with the distance between them. Hence, the magnetic field generated by the transmitter is represented by a magnetic dipole, and the measurements at the Rx are considered to be point measurements of the magnetic field (Farquharson, 2000). Processing routines of data from such small-coil instruments in highly conductive environments cannot be based on existing algorithms because of a nonlinear signal bias introduced by the passive bucking coil (Bx) of the instrument. This coil is used to null the primary field at the receiver location and effectively acts as a second receiver. In low-conductivity environments, the contribution of this second receiver can safely be approximated by a linear factor and the correction is often integrated in the sensor calibration. However, the contribution of the Bx significantly depends on the subsurface conductivity itself (Fitterman, 1998) and needs to be taken into account in the forward model and inversion scheme in highly conductive regimes.

In this study, we applied two modified geophysical inversion algorithms that specifically include the Bx. We used synthetic data and specific instrumental noise characteristics to run inversions for different thickness and conductivity scenarios as a proof of concept for the application to actual field data. We analyzed the sensitivity of the inversion results and their success in resolving sea ice thickness and conductivity simultaneously. Finally, the porosity and its uncertainty is estimated from conductivity retrievals.

We first describe principles behind multifrequency measurements and passive EM bucking, then review two inversion algorithms with the implemented Bx correction for a small-coil EM instrument.

The lower frequencies considered in our study, i.e., $f1=1.53$ and $f2=5.31kHz$⁠, represent measurements in the near-field zone. For instance, for measurements at frequency $f1$⁠, the skin depths in seawater (conductivity of $2.7S/m$⁠) and sea ice of $0.2S/m$ are $δ=7.8$ and $δ=28.8m$⁠, respectively. Both skin depths are significantly larger than the transmitter-receiver separation of 1.66 m. However, measurements at higher frequencies; i.e., $f3=18.33$⁠, $f4=63.03$⁠, and $f5=93.09kHz$ are effectively in the transition zone. For frequency $f5$⁠, $δ=1.0$ and $δ=3.6m$ for seawater and sea ice with a conductivity of $0.2S/m$⁠, respectively. This suggests that in highly conductive sea ice environments, EM soundings by varying frequency are feasible even at transmitter-receiver separations as small as 1.66 m. Moreover, in our case of a moderately conductive layer being underlain by a highly conductive layer, resolution of layer parameters is generally much better than in the opposite case of a highly conductive layer being underlain by a moderately conductive layer (Spies and Frischknecht, 1991).

Current EM induction instruments operate at single to multiple discrete frequencies. A typical EM instrument consists of three coils (Figure 1): Tx, Bx, and Rx. The Tx generates a primary magnetic field inducing a secondary magnetic field in a conductive subsurface. The signal strength of this secondary field depends, among other things, on the distance from the instrument to a body and its conductivity. During a survey, the Rx records the superposition of the primary and secondary magnetic fields minus the bucking field.

The aim of an inversion algorithm is to find a plausible model of the subsurface (usually the simplest), which adequately fits the observed data within the data uncertainty (Farquharson et al., 2003). The process consists of two components: forward modeling and inversion. Forward modeling generates the data of a specific model, whereas inversion automatically changes the model to reduce the misfit between the measured and forward-modeled data (Menke, 1989). We tested two inversion algorithms to resolve synthetic data of homogeneous sea ice for different thicknesses and conductivities. Forward modeling and inversion calculations were performed in 1D. We stitched the 1D inversion results of individual stations, which results in a 2D impression. With the first algorithm (EM1DFM, Farquharson [2000] and Farquharson et al. [2003]), we only inverted for the conductivity within fixed layers of 0.1-m thickness using vertical smoothing between layers. Hence, we call this inversion algorithm smoothness-constrained inversion. Because a vertical conductivity smoothing is not ideal in the case of a sharp sea ice/seawater interface, we applied a second algorithm, the Marquardt-Levenberg inversion (Jupp and Vozoff, 1975; Lines and Treitel, 1984). With only two distinct (instead of many individual) layers and with sea ice thickness as an additional free parameter, this scheme may better reflect reality. In our study, we used the Marquardt-Levenberg inversion algorithm from EM inversion with least intricate algorithms (EMILIA, Kalscheuer et al. [2010, 2012] and Grab [2012]).

The trade-off parameter $β$ (equation 4) determines the balance between the data misfit and the model-structure component and is responsible for fitting the observed data closely, whereas the constructed model should be as simple as possible (Farquharson et al., 2003). The noise is well known. Therefore, we adjusted $β$ with a univariate search until the misfit $ϕd$ matched a predefined target misfit in every iteration $n$⁠. This is the maximum of the number of the calculated responses $Nd=10$ (in-phase and quadrature for the five frequencies), or equal to $0.5ϕdn−1$⁠, to allow only small changes in model structure at each iteration (discrepancy principle, Farquharson, 2000). In case the target misfit cannot be reached, the inversion algorithm searches for the smallest misfit.

To evaluate the compatibility of the new model $mn$ with the field data, the full forward response $dn=F[mn]$ is then used to calculate $ϕd$ (not the approximation in equation 8). The convergence of the inversion algorithm is determined using two convergence criteria according to Gill et al. (1981) and Farquharson (2000).

We applied a model error and resolution analysis to some of our two-layer Marquardt-Levenberg inversions. This method is based on the truncated singular value decomposition (TSVD) of the sensitivity matrix (Kalscheuer and Pedersen, 2007). Because the effective number of model parameters is two (conductivity $σ$ and thickness $t$ of sea ice), the truncation level of the TSVD is set to two in the analyses of these models. Due to the very limited number of model parameters, both layer parameters are perfectly resolved. Hence, in assessing how well the layer parameters are constrained, only the model parameter uncertainties need to be further considered.

To partly account for the nonlinearity of the inversion problem, we computed error estimates for the final inversion models using most-squares inversion (Jackson, 1976; Meju and Hutton, 1992; Meju, 1994). A most-squares inversion is an iterative method that, starting from the least-squares inversion results, determines bounding model values that represent the variability of the model parameters. The goal is to assess the stability of model parameters (Meju and Hutton, 1992). Owing to the logarithmic transformation of model parameters during inversion, the most-squares errors, $1/fMSQ−$ and $fMSQ+$⁠, correspond to parameter ranges $σ/fMSQ−,…,fMSQ+·σ$ and $t/fMSQ−,…,fMSQ+·t$ of layer conductivities and thicknesses, respectively. Hence, a model parameter is well constrained with most-squares errors close to one.

To test the inversion algorithms, we used synthetic data sets for the in-phase and quadrature components, calculated by a 1D forward model using the implementation of Anderson (1979). We calculated two different types of synthetic data sets: (1) neglecting (⁠$dsyn=dRx$⁠) and (2) accounting (⁠$dsyn=dRx−dBx$⁠) for the Bx correction according to equation 3. One data set includes 101 stations, in which the sea ice was increased from 0 to 10 m from one station to the next in increments of 0.1 m. These true thicknesses are hereinafter referred to as $tsyn$⁠. Two exemplary stations are illustrated in Figure 1a, in which a snow layer was neglected. We used typical conductivities $σsyn$ of different sea ice types between 0.01 and $0.2S/m$ (Hunkeler et al., 2015) and a homogeneous half-space conductivity of $2.7S/m$⁠, representing seawater. In-phase and quadrature responses were calculated at frequencies of 1.53, 5.31, 18.33, 63.03, and 93.09 kHz, a setup also used during surveys on sea ice (Hunkeler et al., 2015). For this proof-of-concept of the Bx implementation, we did not account for any noise of the instrument (experiment 1). To simulate realistic data, we added Gaussian noise (corresponding to one standard deviation) to the in-phase and quadrature components based on the free-air response for the respective frequencies (experiments 2–4). The assumed standard deviations between 125 and 204 ppm were calculated by Hunkeler et al. (2015), in which higher frequencies showed higher instrumental noise. The synthetic data from the implementation of Anderson (1979) were compared with the data forward modeled by the inversion algorithms considered here. The data sets computed by all three programs differed only by a few ppm, which is well within the noise level of the instrument.

We applied the two inversion algorithms to synthetic data calculated for stations with variable sea ice thickness and conductivity. We first identify the impact of the Bx on the performance of the EM1DFM inversion (experiment 1). We then investigate whether multilayered, Marquardt-Levenberg inversions are able to reliably reproduce $tsyn$ and $σsyn$ (experiments 2–4).

We used two synthetic data sets, with and without bucking bias, and two versions of EM1DFM, ignoring and considering the bucking bias. We tested the following three scenarios: (1) bucking bias neither in the data nor considered in EM1DFM, (2) bucking bias in the data but not considered in EM1DFM, and (3) bucking bias in the data and considered in EM1DFM. For all calculations, we used $σsyn=0.1S/m$⁠, a seawater conductivity of $2.7S/m$⁠, and $tsyn$ of 0–10 m. For the starting model, we assigned to all stations 1-m-thick sea ice with conductivities of 0.05 and $2.7S/m$ to the underlying seawater. We allowed a maximum of 40 iterations and set the uncertainties in $Wd$ to 1% of the synthetic data.

For the first scenario, the rms error (equation 10) was close to one for most inversions, whereas a rms error $>2$ was found for only three stations (Figure 2a). Although the level of agreement between inversion results and the synthetic data was generally good, an accurate determination of sea ice thickness from these stations was difficult because conductivities varied little between vertically adjacent layers, resulting in a smooth conductivity transition between sea ice and seawater. This transition zone from sea ice to seawater conductivities was small for thinner sea ice (⁠$tsyn<4m$⁠, white steps in Figure 2a), but widened toward thicker sea ice (3 m for $tsyn=10m$⁠). Apart from this variable transition zone, $σsyn$ of $0.1S/m$ (yellow in Figure 2a) and seawater conductivity of $2.7S/m$ (blue in Figure 2a) were reasonably reproduced.

We used for the second scenario the same setup as for the first scenario, and we included the bucking bias in the synthetic data (equation 3; Figure 2b), but not in the inversion algorithm. The algorithm completely failed to reproduce the data for $tsyn<0.5m$⁠, in which high rms errors indicated that the target misfit was often not achieved. For $tsyn$ between 0.5 and 1 m, the algorithm added an additional conductive layer close to the surface. Between 0.5 and 10 m, the transition zone was approximately 1 m too deep compared with $tsyn$⁠. The $σsyn=0.1S/m$ was only rarely reproduced with sea ice conductivity in individual fixed layers being either too high (approximately $0.3S/m$⁠, red in Figure 2b) or too low (close to $0S/m$⁠, white in Figure 2b). For $tsyn>5m$⁠, the algorithm added a conductive layer in the middle of the too-resistive sea ice.

For the third scenario, we used the modified synthetic data from the second scenario together with the modified EM1DFM inversion algorithm (equations 3 and 9). The results (Figure 2c) showed a substantial improvement over the second scenario (Figure 2b), in which the rms error and the resolved sea ice thickness and conductivity were comparable with those from the first scenario. For $tsyn<0.5m$⁠, conductivities similar to the underlying seawater indicate that $σsyn$ could not be resolved.

For all further experiments, we included the Bx in the synthetic data and inversions.

In a second experiment, we used synthetic data with added Gaussian noise and $σsyn=0.01$⁠, 0.05, 0.1, and $0.2S/m$⁠. For each $σsyn$⁠, we calculated the EM response for thicknesses $tsyn$ of 0–10 m. The starting model for each station contained 3 m of $0.05S/m$ (sea ice) and an underlying homogeneous half-space of $2.7S/m$ (seawater). The uncertainty was set to one standard deviation of the noise, and a maximum of 40 iterations was allowed. At all stations, the algorithm stopped before 40 iterations, either because convergence to the target misfit or a minimum was reached, or no suitable model update step was found.

The results of the four scenarios generally showed a good agreement with the synthetic data (Figure 3). The number of stations for which the target misfit could not be achieved were 3, 19, 31, and 29 for $σsyn$ of 0.01, 0.05, 0.1, and $0.2S/m$⁠, respectively. For all inversions, the first two stations with $tsyn$ of 0 and 0.1 m showed high rms errors (⁠$>2$⁠). However, for stations with rms error $<2$⁠, sea ice thickness was reasonably resolved for all four scenarios, but the transition zone between sea ice and seawater conductivity became broader for higher $σsyn$⁠. At station 30 with 3-m-thick sea ice, for example, the zone was for $σsyn=0.01S/m$ only 0.5 m (Figure 3a), whereas for $σsyn=0.2S/m$⁠, it was 1 m wide (Figure 3d). The resulting sea ice conductivity was also more uniform for low $σsyn$ (Figure 3a) compared with the patchy conductivity results from higher $σsyn$ (Figure 3d). However, the sea ice conductivities were clearly distinguishable for the individual scenarios and can be resolved well within $0.05S/m$ using EM1DFM.

For each of the four scenarios, we obtained after the last iteration 10 Jacobian matrices of sensitivities (equation 9) from the five frequencies and both components (in-phase and quadrature). Each Jacobian matrix includes entries at every station in each layer for the respective scenario. To provide a quantified reliability estimate, we normalized all entries of the Jacobian matrices by the respective forward modeled data $d$ obtained after the last iteration. From the calculated normalized Jacobian matrices, we selected for each layer the maximum of the 10 entries. Here, we show for each scenario the resulting 0.5%, 1%, and 2% contour lines, along with the four resulting reliability areas (shaded areas in Figure 3). Small values of normalized Jacobian matrices indicate that a change in the conductivity barely affects the forward modeled data (Farquharson and Oldenburg, 1993); hence, we trust regions with high sensitivities most.

We inverted the same synthetic data used in the second experiment (including Gaussian noise) by the Marquardt-Levenberg algorithm. The data uncertainty was again set to the standard deviation of the noise, and we used a starting model for each station of $0.05S/m$ for the upper 3 m (sea ice) and $2.7S/m$ for the homogeneous half-space below (seawater). The latter was fixed during this inversion because seawater conductivity is usually known and can be determined with tight bounds by sampling through boreholes. In the eastern Weddell Sea, for example, the seawater conductivity ranged from 2.686 to $2.718S/m$ over a time period of roughly three months (Hunkeler et al., 2015) and is usually assumed to be constant during a survey. Tests showed that the Marquardt-Levenberg inversion performed considerably worse when allowing for a free seawater conductivity. Therefore, only sea ice thickness $t$ and sea ice conductivity $σ$ were allowed to vary during a maximum of 300 iterations.

After the maximum of iterations was reached, $t$ was for $tsyn<5m$ well resolved for all four scenarios (Figure 4). The rms error was generally higher for thicker sea ice compared with thinner sea ice in all scenarios. Most stations with a high rms error were found for $σsyn=0.2S/m$⁠, whereas the inversion of approximately half of the stations yielded a $rmserror>1$ (Figure 4d). The lowest rms error and accurately reproduced $t$ (for $tsyn<7m$⁠) were found for the inversions with $σsyn=0.05S/m$⁠, which is also the conductivity defined in the starting model (Figure 4b).

To quantify these findings, we compared the inversion results (and results with $σsyn=0.15S/m$⁠) to $tsyn$ (Figure 5a, scenarios s1–s5). The best resolved sea ice thickness $t$ was found, as stated before, for $σsyn=0.05S/m$ (s2). However, $t$ was generally reasonably reproduced for $tsyn<5m$ in all scenarios.

For comparison to resulting sea ice thickness $t$⁠, we calculated thicknesses with the EMPEX approach. For the latter, we used the 5310 Hz in-phase forward responses, added Gaussian noise, and fitted a single-exponential curve through the data (in-phase response versus sea ice thickness). This curve fit represents the exponential function that is normally calculated from calibration measurements in the field. We used the in-phase component because it is more sensitive to deeper structures, and it is, in contrast to the quadrature component, monotonically decreasing with thickness. The noise of the data was larger than the changes of the exponential function for thick sea ice with low signal responses. Hence, large sea ice thicknesses were poorly constrained. Because the exponential function was not approaching zero, a corresponding sea ice thickness could not always be found for low signal responses and sea ice $>8m$ (Figure 5b).

In general, sea ice thickness $t$ was underestimated with the Marquardt-Levenberg inversion, but overestimated with the calculated exponential function. With both methods, reliable results were found for sea ice $<5m$⁠. For thicker sea ice, the Marquardt-Levenberg inversion depended on $σsyn$⁠, on the starting model and the weighting of the data (more below). The sea ice thickness determination by the exponential approach did not depend on sea ice conductivity because we calculated the exponential fit from the respective in-phase components.

The sea ice conductivity $σ$⁠, however, could be resolved for scenarios s1–s5 (Figure 5a). Figure 5c shows distributions of $σ$⁠, in which individual histograms are clearly separated. For the scenario with the best resolved sea ice thickness $t$ (s2), the narrowest conductivity distribution was found (⁠$0.05S/m$ in Figure 5c). For scenarios s1–s5, the interquartile ranges of $σ$ were within $0.01S/m$ (Figure 5d).

In Table 1, we present the estimated model errors from most-squares inversion of the Marquardt-Levenberg inversion models for stations 11 and 31 (⁠$tsyn$ of 1 and 3 m) for $σsyn=0.05$⁠, 0.1, and $0.2S/m$ (see Figures 4 and 5a). Although the sea ice thicknesses $t$ for station 11 were tightly constrained with $1/fMSQ−=0.99$ and $fMSQ+=1.01$ for all three cases, the sea ice conductivities $σ$ showed more variability with $1/fMSQ−$ and $fMSQ+$ varying from 0.92 to 0.99 and from 1.02 to 1.20, respectively. Similar results were found for thicker sea ice (station 31) with $1/fMSQ−$ between 0.97 and 0.99 and $fMSQ+$ between 1.01 and 1.02 for $t$⁠. The $σ$ showed again more variability with $1/fMSQ−$ between 0.95 and 0.97 and $fMSQ+$ between 1.01 and 1.09. These results were expected because the thickness of a resistive layer is typically better constrained than its conductivity from inductive EM methods (e.g., Pfaffling and Reid, 2009). Similarly, the increase in error with increasing thickness is consistent. The $σ$ is better constrained for higher $σsyn$⁠, whereas the most-squares errors of $σ$ are larger for thinner sea ice. Generally, all errors are very low and all parameters are well constrained because only two model parameters are considered.

To quantify the effect of a priori information, we inverted synthetic data (⁠$σsyn=0.1S/m$ and $tsyn=0−10m$⁠) with different starting models and data weighting.

In a first step, the data set was inverted in a similar manner to experiments 2 and 3. We varied the starting model (Figure 6a) used for each station 3-m-thick sea ice, $σsyn=0.05$ (s1), 0.1 (s2), and $0.2S/m$ (s3) and a homogeneous half-space of $2.7S/m$⁠. The starting model did not have a strong influence on sea ice thickness $t$⁠. The $tsyn<6m$ was reliably reproduced from all conductivities in the starting models (Figure 6a). But even after 300 iterations, the rms error did not reach one (scenarios s1–s3 in Table 2 with rms errors of 1.74–1.87).

In a second step, we used the same starting model as in scenarios s1–s3, but weighted the data misfit according to equation 12. We calculated the weight factors $w$ by using the normalized synthetic data of 5310 Hz in-phase component with added Gaussian noise. Thus, the weight factor $w$ decreased exponentially with increasing station number. Sea ice thickness $t$ (for $tsyn>6m$⁠) was better resolved compared with the results without weighting (Figure 6b), which was also supported by the reduced rms errors of 0.99–1.16 (scenarios s4–s6 in Table 2). For the starting models with sea ice conductivities of $0.05$ and $0.1S/m$⁠, convergence was reached after 104 and 58 iterations, respectively, which was faster than for the other scenarios.

In a third step, we did not weight the data, but we used the sea ice thickness results from the $0.1S/m$ exponential fit from Figure 5b in the starting model. In the cases for which it was not possible during EMPEX processing to find for low in-phase values a corresponding sea ice thickness, we used the thickness from the previous station. This is still apparent in the results, in which layered structures are recognizable for sea ice $>9m$ (Figure 6c). But $t$ again improved for $tsyn>6m$ compared with the case in Figure 6a. The rms errors of 1.04–1.10 for scenarios s7–s9 were smaller than for scenarios s1–s3 (Table 2). Sea ice thicknesses $t$ from scenarios s7 to s9 were similar and therefore did not depend much on sea ice conductivities of the starting model (Figure 6c).

In general, $t$ was reasonably resolved by the Marquardt-Levenberg algorithm for $tsyn<6m$⁠, regardless of weights and starting model. With weight factors and results from single-frequency thickness estimates, $t$ improved for $tsyn>6m$⁠. Although it was not always possible to resolve $t$ up to 10 m, the sea ice conductivities $σ$ were resolved within $±0.01S/m$ for all scenarios (Figure 6d). Outliers of $σ$ rather contained a too-low than a too-high conductivity value.

The cementation factor $m$ depends on the connectivity of the pores and the pore shapes, which are related to the structure of sea ice and the c-axis orientation of individual crystals (Reid et al., 2006). The mobility of ions and the permeability are therefore closely linked to this parameter. With respect to sea ice, the literature suggests $m$ between 1.55 and 2.2 (Thyssen et al., 1974; Morey et al., 1984). Haas et al. (1997) and Worby et al. (1999) use a cementation factor of 1.75 for their calculations. Reid et al. (2006) use values of 1.75 and 2.2, stating that the horizontal conductivity of the sea ice is overestimated with a value of 1.75. Because this parameter is not better constrained for sea ice, and also depends on sea ice types and seasonality, we used a value of 1.75.

With temperature measurements of nine first-year sea ice cores from the winter Weddell Sea (Hunkeler et al., 2015), we calculated a mean brine conductivity of $4.69±0.91S/m$ according to Stogryn and Desargant (1985). At the same locations, Hunkeler et al. (2015) obtain from EM measurements at 63.03 and 93.09 kHz an average bulk conductivity of $0.06S/m$ from seven experiments above thin sea ice (⁠$<1m$⁠). Assuming now $σb=4.69±0.91S/m$⁠, $m=1.75$⁠, and $σ=0.06±0.01S/m$⁠, we obtain a porosity of $ϕ=8.3±1.2%$ using Gaussian law of error propagation. From the nine sea ice cores, we calculated according to Cox and Weeks (1983) and Leppäranta and Manninen (1988) average brine volumes (porosities) of $ϕ=9.5±5.0%$⁠. A reasonable sea ice porosity can be estimated from bulk conductivity obtained by EM, according to the agreement between porosities calculated from conductivities and direct measurements in sea ice cores.

We have shown that 1D inversions of EM data from passively bucked small coil systems must include the nonlinear contribution of the Bx when working close to a highly conductive environment. In our case of sea ice thickness and conductivity retrieval, the negligence of this contribution would lead to wrong results for sea ice thinner than 2 m (Figure 2b), which is typical of level sea ice thickness in the Arctic and Antarctic landfast sea ice. We calculated synthetic data based on a typical commercial instrument used in sea ice research (Figure 1). The implementation of the bucking bias, however, will need to be made for all frequency-domain EM sensors that use passive bucking in highly conductive environments (Fitterman, 1998). With the effect included, we were able to resolve conductivity and ultimately porosity with uncertainty values that motivate further studies and use of 1D inversions for field data. This and earlier case studies have shown that, especially, a high-frequency range is required for the retrieval of sea ice conductivity.

Our aim was to jointly resolve sea ice thickness $t$ and conductivity $σ$ for a single layer with a thickness range of $tsyn$ of 0–10 m and conductivities $σsyn$ of $0.01−0.2S/m$⁠. We have tested two inversion methods, a smoothness-constrained (EM1DFM) and a Marquardt-Levenberg (EMILIA) inversion. In each scenario, we inverted synthetic data with constant conductivity $σsyn$ and the full range of physically feasible sea ice thicknesses $tsyn$⁠, in which the upper range includes extreme cases, such as decadal sea ice. During each iteration, one single $β$ (trade-off parameter [EM1DFM] or damping factor [EMILIA]) was used for all stations simultaneously. The alternative would be separate inversions of all stations with individual $β$⁠. This approach may lead to a faster convergence with fewer iterations. However, we chose the single $β$ approach under the assumption that individual stations are not independent to reduce the potential for a wrong minimum of individual inversions.

EM1DFM inverts for the conductivity of multiple fixed layers; therefore, the resolution of thickness depends on how the vertical conductivity gradient is evaluated (Figure 3). The conductivity transition between sea ice and seawater becomes broader for thicker sea ice, and it also depends on $σsyn$⁠, thus subsequently increasing the uncertainty of the thickness estimate. The Marquardt-Levenberg algorithm instead uses layer thickness and layer conductivity as free parameters, therefore removing the need for retrieving layer interfaces from conductivity gradients (Figure 4). Although the description of sea ice as a layer of defined thickness with a conductivity discontinuity at the ice-water interface is a more realistic description of undeformed sea ice, internal structures such as from sea ice deformation or complex layering potentially can be better described by EM1DFM. However, the uncertainty of the interface location will be increased by the process of vertical conductivity smoothing by EM1DFM. The choice of the inversion scheme therefore depends on the user and the specific application.

The best sea ice thickness results from the Marquardt-Levenberg inversions were found using data weights or tuning the starting model using single-frequency processing results (Figure 6). Higher weights for the weaker signals of thicker sea ice do increase the convergence rate and precision (Table 2). For real field data, we need to be careful when weighting low signal responses because signal noise and instrument attitude errors may have an increased effect. In any case, with both inversion methods, $tsyn<5m$ were reasonably resolved. The conductivities of the starting model seemed to have little effect on the inversion results (⁠$tsyn<5m$⁠). The constant 3-m-thick starting probably favored the better retrieval of thin sea ice.

Although the sea ice thickness is better constrained than the sea ice conductivity based on the most-squares inversion (Table 1), the conductivity is determined by the Marquardt-Levenberg inversion with the five used frequencies within $±0.01S/m$ no matter how well the thickness was resolved, independent of the starting model, the weight factors, and $σsyn$ (Figures 5d and 6d). Likewise, sea ice conductivities can be reasonably well resolved for most of the fixed layers of EM1DFM inversions, provided that the rms error $<2$ (Figure 3).

We have simplified sea ice as one homogeneous layer of prescribed thicknesses and conductivities, not accounting for any anisotropy due to higher conductivities in vertical brine channels (Reid et al., 2006). In nature, the simplification of level sea ice most likely applies to thermodynamically grown sea ice, which grows to a maximum of 2.5 m (Thomas and Dieckmann, 2009) or to the Antarctic subice platelet layer with multiple layers of different conductivities (Hunkeler et al., 2015). The sea ice/ocean interface or the interface between sea ice and the platelet layer is characterized by discontinuities in the vertical conductivity structure and by far smaller variability in the horizontal direction. To resolve such structures, we will put effort into improving laterally constrained inversion schemes (Auken et al., 2005). But even such pseudo-2D inversions cannot account for sea ice deformation features, in which conductivities and subice topographies vary additionally on a subfootprint scale. To account for complex topographies of deformed sea ice, we need to include the bucking bias in a full 2D inversion, in which forward-modeling responses and sensitivities are calculated in 2D.

The conversion of conductivity into porosity requires the knowledge of two parameters in Archie’s law: (1) the cementation factor and (2) the brine conductivity. First, we assume a cementation factor that is based on the results from previous sea ice field experiments. Second, we used brine conductivities found in Antarctic winter sea ice cores.

For sea ice throughout the freezing season, this parameter is of second-order significance because it is not assumed to vary much. In summer, the brine contains a mixture of saline brine and freshwater due to the surface and internal melting, yielding a higher porosity than winter or spring sea ice (Vancoppenolle et al., 2007; Tison et al., 2008). In such cases, we can still make assumptions about the brine conductivity, but because we deal with higher uncertainties, the effect of porosity and brine conductivity on the bulk sea ice conductivity might not be fully separated in summer cases.

However, the uncertainty range of the sea ice conductivity in this idealized case is promising for retrieving porosity from actual field data. The improvement of 1D sea ice inversions presented in this study therefore has the potential to improve biogeochemical sea ice studies (e.g., Arrigo et al., 1997) and estimations of sea ice density for satellite altimetry retrieval of sea ice thickness (Ricker et al., 2014) by large-scale mapping of sea ice porosity.

Sea ice conductivity is a proxy parameter that can be used to derive porosity, albeit systematic measurements on scales larger than point measurements have been rare mostly due to a lack of suitable methodology. With the nondestructive method of multifrequency EM sounding, this study aims to lay the foundation for an efficient joint retrieval of sea ice thickness and porosity. We varied parameters of a synthetic data set (thickness and conductivity) to assess the sensitivity of different parameter combinations in two geophysical inversion algorithms. Furthermore, we tested the influence of the starting model and data weighting. We calculated the synthetic data for a small-coil system with frequencies relevant for sea ice research in the range of 1–100 kHz. The existing inversion algorithms were modified for a nonlinear correction of signal bias that is caused by passive bucking of the primary EM field close to highly conductive environments. Only with the integration of this sensor-specific correction were we able to determine sea ice conductivity with a resolution of $0.01S/m$ with 1D inversion methods. This resolution is sufficient to initiate field trials to estimate sea ice porosity of different sea ice types, the thickness and porosity in multiyear cases, flooded snow layers, or subice platelet layers near Antarctic ice shelves. The detailed knowledge of porosity facilitates biogeochemical studies of sea ice and sea ice densities estimated from porosities yield freeboard validations of satellite measurements. Although our approach improves sea ice surveying with EM induction sounding methods, it suffers the same limitations caused by the use of 1D forward models for essentially 3D targets with a significant subfootprint-scale variability. Because real-world sea ice surveys are mostly conducted along transects, we suggest modifying 2D inversions methods to resolve thickness and internal conductivity changes at subfootprint scales. Such a method would enable further studies of sea ice deformation, relevant for the improved estimation of sea ice volume.

Computing resources were provided by SNIC through the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX, Uppsala, Sweden) under project snic2014-1-243 and by Eidgenössische Technische Hochschule (ETH) Zurich, Switzerland. J. Hermansson at UPPMAX and H. Horstmayer at ETH are acknowledged for assistance concerning technical and implementational aspects in making the code run. This work was partly funded by the POLMAR graduate school and the Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung. We are very grateful to the reviewers whose comments and suggestions improved the clarity of the manuscript.