## Abstract

An analytic continuation method for obtaining rigorous bounds on the effective complex permittivity ** ε*** of polycrystalline composite materials is developed. It is assumed that the composite consists of many identical anisotropic crystals, each with a unique orientation. The key step in obtaining the bounds involves deriving an integral representation for

***, which separates parameter information from geometrical information. Forward bounds are then found using knowledge of the single crystal permittivity tensor and mean crystal orientation. Inverse bounds are also developed, which recover information about the mean crystal orientation from**

*ε****. We apply the polycrystalline bounds to sea ice, a critical component of the climate system. Different ice types, which result from different growth conditions, have different crystal orientation and size statistics. These characteristics significantly influence the fluid transport properties of sea ice, which control many geophysical and biogeochemical processes important to the climate and polar ecosystems. Using a two-scale homogenization scheme, where the single crystal tensor is numerically computed, forward bounds for sea ice are obtained and are in excellent agreement with columnar sea ice data. Furthermore, the inverse bounds are also applied to sea ice, helping to lay the groundwork for determining ice type using remote sensing techniques.**

*ε*## 1. Introduction

A polycrystalline composite material consists of many single crystals that can vary in shape, size and orientation. A broad range of manufactured and naturally occurring materials are polycrystalline, including metals, ceramics, rocks, glacial ice and sea ice. Here, we consider the electromagnetic behaviour of polycrystalline media when the wavelength is much larger than the scale of the underlying microstructure of the composite. When in this regime, the quasi-static approximation is valid, and the electric and displacement fields can be viewed as time-independent fields. Then the polycrystalline composite can be characterized electromagnetically via the effective complex permittivity tensor ** ε***.

The macroscopic permittivity or dielectric tensor ** ε*** of a polycrystalline composite depends directly on its microstructural properties, such as the complex permittivity tensor of the individual crystals and their microstructural geometry, i.e. how the crystals are oriented. Owing to the complicated nature of the microstructure, explicitly calculating

*** is highly non-trivial, and can generally only be accomplished if the exact microstructure is known and with the assistance of very powerful numerical computations. Therefore, using partial microstructural information that may be available to estimate or bound**

*ε**** is a very practical and useful approach.**

*ε*There has been extensive work in the past on estimating and bounding ** ε*** for composite materials. The books by Cherkaev [1] and Milton [2] thoroughly discuss much of this work. In particular,

*** has been intensively studied for two phase composites. Rigorous bounds were first obtained in the early 1980s using the analytic continuation method, where the effective parameter is treated as an analytic function of the ratio of the component parameters [3–6]. These bounds assume that the complex permittivity of each component is known and that there is some partial information available about the microstructure. The most general bounds assume only knowledge of the relative volume fractions of each material, resulting in the complex versions of the classical arithmetic and harmonic mean bounds for a two-component material. Tighter bounds can be found when more geometrical information is available, such as knowing that the microstructure is isotropic [7], or that the composite has a**

*ε**matrix–particle*structure [8,9], etc.

Additionally, the electromagnetic response of a composite material can be used to help determine microstructural properties when approached as an inverse problem. That is, given information on ** ε***, different microstructural details can be resolved, such as the relative volume fractions of each component of the material. This has also been extensively investigated for a two-component composite [10–21], and computational approaches [15,16,18–20] as well as analytic inverse bounds for geometric parameters [12–14] have been developed.

For some composite structures, it is more appropriate to assume that the material consists of many identical anisotropic pieces that are oriented in different directions. This is the case for a polycrystalline composite. Polycrystalline materials have been studied for decades [1,2,22–30], and in particular, there has been a significant amount of work done on bounding the (real) effective conductivity. The books by Cherkaev [1] and Milton [2] discuss the majority of this work.

Here we develop an analytic continuation method for obtaining complex bounds on ** ε*** for a three-dimensional (transversely isotropic or uniaxial) polycrystalline composite material. The key step in obtaining the bounds involves deriving an integral representation [3,6] for the effective complex permittivity tensor. This approach was first employed by Bergman [3], and the bounds for the diagonal elements of the complex permittivity tensor

*** for a two-component composite were derived in [3–6,31] using this approach. An important feature of the integral representation is that it separates parameter information from geometrical information. By making an assumption about the complex permittivity tensor of each individual crystal and assuming some knowledge about the mean single crystal orientation, we obtain first-order polycrystalline bounds on**

*ε**** for the entire polycrystal. If we further assume that the polycrystalline material has the ‘polycrystalline Hashin–Shtrikman condition’ [2], which is essentially geometric isotropy, second-order forward bounds are constructed. Further, we use an inverse analytic bounds method [12–14] and derive inverse bounds for the mean orientation of crystals in the polycrystalline composite. Thus, knowing**

*ε**** and the complex permittivity tensor of an individual crystal, we bound the mean single crystal orientation.**

*ε*As a demonstration of the complex polycrystalline bounds, we compare them to sea ice data. Investigating the electromagnetic behaviour of sea ice is not only interesting from a composite material point of view, but also because of the valuable information that can be recovered using remote sensing techniques, such as sea ice thickness and fluid transport properties. Sea ice covers between 7 and 10% of the Earth's ocean surface and is both an indicator and agent of climate change [32–34]. Since the 1980s, there has been a steady decline in Arctic summer sea ice extent, with a much more rapid decline over the past decade [35]. During the winter months in the Arctic and Antarctic, the extensive sea ice packs serve as the boundary layer which mediates the exchange of heat, moisture and momentum between the atmosphere and ocean [33,36]. The vast expanses of sea ice also serve as a habitat for rich microbial communities living in the brine microstructure of porous sea ice [33,37,38]. In turn, these microbial communities are primary providers for the complex food webs in the polar oceans.

Owing to the global nature of monitoring the Earth's sea ice packs, large-scale information is usually obtained via remote sensing from platforms on satellites, aircraft and ships [17,39–42]. One of the grand challenges of sea ice remote sensing is to accurately recover the thickness distribution of the pack. Assessing the impact of climate change on the polar regions involves monitoring not only the ice extent, but the ice volume, which requires knowledge of ice thickness. Recently, there has been increasing interest in using low-frequency electromagnetic induction devices to estimate sea ice thickness [43]. In addition to assessing ice thickness, remotely monitoring the fluid transport properties of sea ice is of increasing interest because of the broad range of geophysical and biological processes it mediates in the polar marine environment. For example, the evolution of melt ponds and summer ice albedo is constrained by drainage through porous sea ice [44], where ice-albedo feedback is believed to play a key role in the decline of summer Arctic sea ice [35]. Fluid flow also facilitates snow-ice formation [45], the evolution of the salt budget [33], convection-enhanced thermal transport [46], CO_{2} exchange [47] and biomass build-up sustained by nutrient fluxes [33,37]. There is evidence [48] that the polycrystalline structure of the sea ice, such as granular versus columnar, can dramatically affect its fluid transport properties. Thus, determining ice type using remote sensing techniques may be a particularly useful application.

There has been considerable work in the past on estimating and bounding ** ε*** for sea ice, particularly in the microwave region [9,14,41,49–57]. The rigorous two-component bounds mentioned above have successfully been used to bound

*** for sea ice [9,14,54–57]. These bounds assume that sea ice is a two-component material, consisting of a pure ice and brine phase. The forward bounds recover information on**

*ε**** using information about the microstructure, such as brine volume fractions or porosity**

*ε**ϕ*(and sometimes further assuming statistical isotropy), while the inverse bounds attempt to recover

*ϕ*from

***.**

*ε*Here we apply the first-order forward polycrystalline bounds to sea ice. We see a dramatic improvement over the classic two-component bounds, because these new bounds include additional information about single crystal orientations. Here, we use the dataset presented in [52] to compare the polycrystalline bounds to sea ice. This dataset is the same one used in [9,14,54], thus helping provide some continuity between different types of bounds. In addition to providing ** ε*** and

*ϕ*measurements, the dataset provides detailed crystallographic data, which will be critical when applying the bounds. Notationally, we will reserve

*R*

_{1}and

*R*

_{2}to indicate the previously reported two-component forward bounds and use

*R*

_{3}and

*R*

_{4}to describe the new polycrystalline forward bounds. The single crystal complex permittivity tensor is obtained by numerical simulation using X-ray CT data on sea ice, along with known brine volume fractions and ice and brine permittivities. A governing single crystal complex permittivity tensor, that is applicable to an entire sea ice column, is then obtained using several generalizing assumptions and approximations. Further, the inverse method that we develop is applied to sea ice and we obtain bounds on the mean single crystal orientation. Columnar and granular microstructures (figure 1) have different mean single crystal orientations [58], thus this inverse approach helps lay the groundwork for determining ice type when using remote sensing techniques.

## 2. Forward bounds on the effective complex permittivity of a polycrystalline material

Consider the constitutive relation **D**(**x**,*ω*)=** ε**(

**x**,

*ω*)

**E**(

**x**,

*ω*), where

**D**(

**x**,

*ω*) and

**E**(

**x**,

*ω*) are stationary random displacement and electric fields and

**(**

*ε***x**,

*ω*) is the permittivity tensor of some medium. Here

*ω*∈

*Ω*, where

*d*is the spatial dimension and

*Ω*is the set of all realizations of the random medium. Let us consider a polycrystalline material, where each crystal has the same complex permittivity tensor but with different orientation. Thus, we let

**(**

*ε***x**,

*ω*)=

**B**

^{−1}(

**x**,

*ω*)

*ε*_{d}

**B**(

**x**,

*ω*), where

**B**(

**x**,

*ω*) is a rotation matrix describing the orientation of a crystal at location

**x**and realization

**, and**

*ω*

*ε*_{d}is the same permittivity tensor for each crystal and can be written

Here, we are assuming that each crystal has the same permittivity value *ε*_{2} in both horizontal directions (transversely isotropic or uniaxial) with anisotropy occurring in the vertical direction with permittivity value *ε*_{1}. It is assumed that *ε*_{1} and *ε*_{2} can take complex values. Then making the assumptions that we are in the quasi-static regime and there is no free charge, we can write **∇**×**E**(**x**,*ω*)=0 and **∇**⋅**D**(**x**,*ω*)=0. Now, letting 〈⋅〉 represent an ensemble average over *Ω* or a spatial average over all of **E**(**x**,*ω*)〉=**e**_{k}, where **e**_{k} is a unit vector in the *k*th direction for some *k*=1,…,*d*. For notational simplicity, we write **E**(**x**,*ω*)=**E** and **D**(**x**,*ω*)=**D**. The effective complex permittivity tensor is then defined via

From this, we can write ** ε***=[

***]**

*ε*_{kk}. This allows us to strictly examine the

*kk*th component of the effective permittivity tensor and simplifies the notation. Thus, we can rewrite the equation as

***(λ**

*ε**ε*

_{1},λ

*ε*

_{2})=λ

***(**

*ε**ε*

_{1},

*ε*

_{2}),

*** depends only on the ratio**

*ε**h*=

*ε*

_{1}/

*ε*

_{2}and we define

*m*(

*h*)=

*ε**/

*ε*

_{2}. Therefore, we have the equation

**C**=

**e**

_{1}(

**e**

_{1})

^{T},

**I**is a 3×3 identity matrix, and

**e**

_{1}is a unit vector in the first direction [2,29]. To simplify the notation, we define

**R**=

**B**

^{−1}

**C**

**B**, and can then write

The two main properties of *m*(*h*) are that it is analytic off *h*-plane, and that it maps the upper half plane to the upper half plane [3,6], so that it is an example of a Herglotz or Stieltjes function [59]. The key step for obtaining forward bounds is to use an analytic continuation method which involves obtaining an integral representation for ** ε***. If we let

*s*=1/(1−

*h*), then we can define

**E**, which will allow us to find an integral representation for

*F*(

*s*).

To find the resolvent representation of **E**, first examine **∇**⋅**D**=0, which implies that **∇**⋅*ε***E**=0. Then, let **G** be a vector representing the mean fluctuations in the electric field and call **E**=**e**_{k}+**G**. Expand ** ε** and

**E**using the previous definitions and formulate the equation

**I**is again a 3×3 identity matrix. After multiple algebraic manipulations, one can obtain the equation

**as the operator**

*Γ***∇**(−Δ)

^{−1}

**∇**⋅, which projects fields onto a subspace of curl-free mean-zero fields

**in [62]. The same operator is also used in [16,54,56,57]. Using it, we find the resolvent representation for**

*Γ***E**to be

*F*(

*s*) using the following equation:

*F*(

*s*) which takes the particularly nice form

*μ*on [0,1] is the spectral measure of the self-adjoint operator

*Γ***R**.

*F*(

*s*) is also analytic off [0,1] in the

*s*-plane, which is the only restriction for this integral representation. All of the geometrical information is now contained inside of

*μ*and all of the parameter information is contained in

*s*, including the electromagnetic wave frequency. Expanding

*F*(

*s*), we find

*μ*via its moments

*k*th direction.

Bounds on *ε**, or *F*(*s*), are obtained by fixing *s* in (2.7), varying over admissible measures *μ* (or admissible geometries), such as those that satisfy only *F*(*s*) in the complex plane [6]. The bound *R*_{3} assumes only that the mean crystal orientation *μ*, the extreme values of *F* are attained by extreme points of *ε** lie inside the region *R*_{3} bounded by circular arcs, one of which is parametrized in the *F*-plane by
*E*(*s*)=1−*ε*_{1}/*ε**, which is a Herglotz function like *F*(*s*), analytic off [0,1]. Then in the *E*-plane, we can parametrize the other circular boundary of *R*_{3} by
*ε**-plane, *R*_{3} takes the following form for the ‘lower’ and ‘upper’ bounds, which are still circular arcs, *ε**_{l} and *ε**_{u}, respectively.

When *ε*_{1} and *ε*_{2} are real and positive, the bounds collapse to the interval
*k*.

To obtain second-order complex bounds further assumptions need to be made. For instance, if the polycrystalline composite is assumed to have the ‘polycrystalline Hashin–Shtrikman condition’ [2] or essentially geometric isotropy, then *d* is the dimension of the polycrystalline composite. In two dimensions, we define a polycrystalline material to be geometrically isotropic if for every crystal in the polycrystalline composite with orientation off the vertical direction described by the normalized vector 〈*x*,*y*〉, there exist three other crystals that have orientations 〈*x*,−*y*〉, 〈*y*,*x*〉 and 〈*y*,−*x*〉. A similar definition can be made for three dimensions, where groups of 24 crystals are needed instead of groups of four. (Note: several special examples in two dimensions where only groups of two are required include polycrystalline materials where all the single crystals are vertically and horizontally aligned or all the single crystals have an orientation angle of ±*π*/4 radians or ±45° off the vertical axis. In a similar manner, if all the single crystals are vertically or horizontally aligned in three dimensions, only groups of three crystals are required.)

Here we show the derivation for the value of ** Γ**=

**∇**(−Δ)

^{−1}

**∇**⋅ and define (−Δ)

^{−1}in terms of a Green's function so that

*g*(

*x*)=−

*δ*

_{y}(

*x*) and

**R**takes the form

*θ*is the angle of orientation off the vertical axis. Define

Again for simplicity, let us consider the terms separately. The intent here is to find crystal orientation combinations so that the Laplacian operator is recovered, and together with *g*(*x*,*y*) we obtain a delta function inside the integral, or for the term to become zero. Examining the first term, analysing the first direction (*k*=1), and using the notation *D*_{x}=d/d*x*, we see that
*a*^{2}+*b*^{2})*D*_{xx}−2(*ab*+*bc*)*D*_{xy}+(*b*^{2}+*c*^{2})*D*_{yy}] to become the Laplacian operator, *a*^{2}=*c*^{2} and *b*=0. This is equivalent to every crystal in the polycrystalline material having either perfect vertical or horizontal rotations with an equal amount of crystals in the vertical and horizontal directions. Therefore,

This same line of reasoning can be expanded into a geometrically isotropic polycrystalline material. That is, let us now consider ‘groups’ of four single crystals under the statistical average. If four crystals are examined at once, the condition for recovering the Laplacian operator becomes
*q* is some constant. Thus, *a*_{1}*b*_{1}+*a*_{2}*b*_{2}+*a*_{4}3*b*_{3}+*a*_{4}*b*_{4}+*b*_{1}*c*_{1}+*b*_{2}*c*_{2}+*b*_{3}*c*_{3}+*b*_{4}*c*_{4}=0. These conditions are satisfied provided that the orientation off the vertical directions for the four crystals, in terms of normalized vectors are 〈*x*,*y*〉, 〈*x*,−*y*〉, 〈*y*,−*x*〉 and 〈*y*,*x*〉. Conveniently, for these four specific crystals, *q* takes the value of *q*=1/*d*. Therefore, the same result holds for the statistical average as before and we see that

Now consider the second term from equation (2.15), which is *k*=1,2), then over the statistical average
*a*_{1}+*a*_{2}+*a*_{3}+*a*_{4}=*c*_{1}+*c*_{2}+*c*_{3}+*c*_{4} and *b*_{1}+*b*_{2}+*b*_{3}+*b*_{4}=0. Therefore, the second term also introduces a Laplacian operator and we see that

Finally, consider the third term from equation (2.15), which is

Putting all together, we see that *μ*_{1} in equation (2.15), under the assumption of geometrical isotropy, satisfies
*μ*_{1} is analogous to the

Thus, if the polycrystalline material is further assumed to have the Hashin–Shtrikman condition, then *F*(*s*) is known to second order, with *μ*_{0}=1/*d* and *μ*_{1}=(*d*−1)/*d*^{3}, so that

A convenient transform *F*_{1}(*s*) is an upper half plane function analytic off [0,1] and has the representation

Under the additional assumption of geometric isotropy *F*_{1}(*s*) is known only to first order, where *F*_{1}(*s*)=(*d*−1)/(*ds*)+⋯ , and *F*_{1}(*s*) lie in the circular arc (*d*−1)/[(*d*)(*s*−*z*)], *F*-plane, we can parametrize one boundary of *R*_{4} by

Similarly, to display the other arc, we use the auxiliary function *E*(*s*)=1−*ε*_{1}/*ε**=(1−*sF*(*s*))/(*s*(1−*F*(*s*))), and find that *E*(*s*)=(*d*−1)/(*ds*)+(*d*^{2}−2*d*−1)/(*d*^{3}*s*^{2}). Again, using a similar method as with *F*(*s*), an arc can be found in the *E*-plane, and we can parametrize the other circular boundary of *R*_{4} by

When *ε*_{1} and *ε*_{2} are real and positive, the bounds collapse to the interval
*s*=*ε*_{2}/(*ε*_{2}−*ε*_{1}) and *ε*_{2}≤*ε*_{1}. These are exactly the Hashin–Shtrikman bounds for an isotropic polycrystalline composite [22]. Further, if we evaluate the two-dimensional second-order complex forward bounds for a two-component material [4–6], where each material has a volume fraction of 50%, we see that they are in agreement with the two-dimensional second-order complex forward polycrystalline bounds. The three-dimensional bounds are also in agreement.

For the purpose of comparing these bounds to previously established ones on the (real) effective permittivity [1,2,23–27], let us further consider that the polycrystal is isotropic in the sense that *et al.* [25] is achieved with the *sphere assemblage model* conjectured by Schulgasser [23,24]. The reason being that the conductivity in each direction is simultaneously minimized in the equation ** ε***=(1/3) tr(

***)≥**

*ε**ε*

_{s}[25], where

*ε*

_{s}is the permittivity of the

*sphere assemblage model*and

*** is the full permittivity tensor. Therefore, a minimum is found and achieved with the**

*ε**sphere assemblage model*because all directions have the same minimum permittivity. The lower bound we find here is only for

## 3. Inverse bounds for structural parameters

The objective of inverse bounds is to use data from the electromagnetic response of a polycrystalline material to recover information about its structural parameters. In previous work [10–14,17], this is typically done to recover information about the volume fractions of the two constituents of a composite material. Here, we will show how to recover information about the mean crystal orientation *ε**, *ε** touches one boundary of the region *R*_{3} described in the previous section and is then decreased until the value touches the other boundary. This procedure gives an analytic estimate (the first-order inverse bounds [12–14]) of the range of values of the mean crystal orientation *f* is the known value of *F*(*s*) and *g* is the known value of *G*(*t*)=1−*ε**/*ε*_{1} with *t*=1−*s*.

The objective of the second-order inverse bounds would be to obtain a better estimate for the mean orientation of crystals in the *k*th direction. However, as demonstrated in the second-order forward bounds, the mean orientation must be the same in all directions. Thus, we already know that,

## 4. Comparison of the bounds to sea ice data

Here the polycrystalline bounds derived in the current work are applied to sea ice composite, and the results are compared with the measured effective permittivity of sea ice in [52]. This dataset is obtained from primarily columnar sea ice; it was previously used to compare effective permittivity of sea ice with the bounds for the effective property of a two-component material [14,54], of a statistically isotropic two-component composite [14,54] and of a two-component matrix–particle material [9,21]. Applying the polycrystalline bounds to the same set of data allows for comparison between different types of bounds and a deeper understanding of the physical relationship the polycrystalline bounds provide in the case of sea ice data. To calculate the forward polycrystalline bounds, the method uses information on the complex permittivity tensor of the identical single crystals (i.e. *ε*_{1} and *ε*_{2}) and the crystal orientation statistics (i.e.

The single crystal complex permittivity tensor for sea ice is obtained by evaluating X-ray CT data with known ice and brine permittivities and brine volume fractions *ϕ* using Comsol 3.5a. We examine 222 single crystals at a frequency of 4.75 GHz and at a temperature of −6°*C*, where brine has a permittivity value of 51.07+45.160*i* [64] and pure ice has a permittivity of 3.15+0.002i. In reality, ice is a three-component composite consisting of pure ice with air bubbles and pockets of brine. As was demonstrated in [14], it is important to model the effect of the air phase in sea ice when calculating the complex permittivity. Thus, to account for the air phase, a Maxwell-Garnett mixing formula is used to calculate the permittivity of the ice with air bubbles as was done in [14]. The air phase is disconnected, has a small volume fraction and is contained entirely within the ice phase, making the Maxwell-Garnett mixing formula an extremely accurate approximation for the permittivity of the combined air–ice phase. Therefore, the permittivity used for the air–ice phase is 3.07+0.002*i*. Different single crystal microstructures were calculated at different values of brine volume fractions *ϕ* and a dataset of single crystal complex permittivity tensors was generated (table 1).

Different sea ice single crystal geometric configurations can have significantly different permittivity tensors for the same brine volume fraction *ϕ* value as is shown in table 1. The differences observed for single crystals with the same brine volume fractions are largely due to the different geometric configurations of the brine phase within the sample. The samples were taken at a temperature of −6°*C*, which is close to the percolation threshold, and the effect of the changes in the microgeometry could be very pronounced. This is in some sense a finite-size effect but is expected due to the different samples having finite sizes. In particular, the permittivity in the vertical direction can dramatically change depending on the brine connectedness in the vertical direction. Furthermore, the two horizontal components tend to have slightly different permittivity values. The polycrystalline bounds are derived under the assumption that the composite material is transversely isotropic or uniaxial, i.e. it is composed of many identical crystals with the same permittivity in two (horizontal) directions and a different permittivity in the other (vertical) direction. As is quickly observed in table 1 actual sea ice does not exactly satisfy these assumptions, namely ice crystals vary and the horizontal permittivities differ from each other. Although the fraction of brine *ϕ* within a single crystal displays little variance from the top of the crystal to the bottom of the crystal, the fraction of brine *ϕ* can dramatically change across an entire sea ice column, thus substantially changing the single crystal permittivities at different depths. For example, the very bottom layer of a sea ice column can have a brine fraction *ϕ* almost twice as large as the average, which is typical in classic columnar sea ice [65,66]. Further, as displayed in table 1, the change of permittivity of the vertical component is not linear with respect to *ϕ*, and averaging permittivities over a small range of *ϕ* values to obtain a single crystal permittivity tensor for the entire ice column will not accurately represent the physics of the ice.

To account for these differences between the assumptions under the polycrystalline bounds and actual sea ice data, we will ‘idealize’ the sea ice data, so that the polycrystalline bounds may still be applied. Inherently, this idealization changes the objective from finding an exact forward bound for a specific configuration with identical single crystals to finding a more general forward bound that can be applied to a large class of sea ice, such as all columnar sea ice within a certain brine volume fraction *ϕ* range. As part of this idealization, a governing single crystal permittivity must be found that reflects the overall effect of the different single crystal permittivity tensors that span the entire ice column. Without detailed information on the entire ice column, getting an exact permittivity for the governing single crystal is incredibly difficult if not impossible. However, a plausible approximation for the governing single crystal can be found by averaging different single crystal permittivities. Thus, to obtain the governing single crystal permittivity tensor for a value of *ϕ* corresponding to the entire ice column, we averaged the permittivities using arithmetic means in the vertical direction and both of the horizontal directions over a range of *ϕ*, accounting for a typical variation in *ϕ* across an entire ice column. For example, for an average value of *ϕ*=3.5% we used the following range of (averaged) *ϕ* values: 0.025, 0.025. 0.025, 0.03 and 0.07. For an average value of *ϕ*=4%, we used the following range of (averaged) *ϕ* values: 0.03, 0.03, 0.03, 0.03, 0.04 and 0.08. Therefore, the governing single crystal permittivity tensor for an entire column of sea ice with overall values of *ϕ*=3.5% and *ϕ*=4% are: *ε*_{1}=3.74+0.62*i* (vertical), *ε*_{2}=3.46+0.08*i* (horizontal) and *ε*_{1}=4.11+0.67*i* (vertical), *ε*_{2}=3.52+0.10*i* (horizontal), respectively.

The polycrystalline forward bounds also incorporate information on the single crystal orientation statistics *c*-axis distribution statistics of sea ice as a function of depth for ice grown in a region without a preferred current direction (thus, transversely isotropic). The big picture is that granular ice (typically found in the top layer of a column of sea ice [58]) has an essentially random uniform distribution across all angles, whereas columnar ice has a strongly preferred vertical orientation. The effective complex permittivity dataset from Arcone *et al*. [52] is largely columnar and cross-sectional slices show that the sea ice is transversely isotropic. Examining the orientation statistics [58], it is very reasonable to conclude that the average crystal orientation

Owing to the necessary idealizations describe above, the objective of the forward bounds is to capture all possible effective complex permittivity variations that can occur in columnar ice. The electromagnetic wave propagating through the sea ice column in [52] is orthogonal to the horizontal plane, thus the electric field is in the horizontal plane. The sea ice structure in the horizontal plane is isotropic, therefore, we examine one of the horizontal directions (depending on the direction of the applied electric field). This allows us to reduce the three-dimensional problem to a two-dimensional one and to use a two-dimensional rotation matrix, which gives the same result as a three-dimensional rotation matrix in the case when the applied electric field is in one of the horizontal directions. Therefore, if the electric field is in the *k*=2 horizontal direction, for an average deviation off the vertical axis between 0° and 30° (i.e. columnar ice), the crystal orientation

The first-order polycrystalline forward bounds can then be applied to the data and the largest area of overlap between the bounds is assumed to be the region where the data must lie. It is possible that the forward bounds overestimate the region because of this technique (namely, we assume the mean orientation is between 0° and 30° off the vertical axis). However, each region in these bounds could still be found by slightly adjusting the single crystal permittivity tensor (which can have some variability) for a specific orientation. Further justification of this approach is the large (and possibly unknown) variability in the general columnar crystal orientation statistics from sample to sample. The bounds are general enough to accurately predict the permittivity of primarily columnar sea ice without having to know specific orientation statistics or a specific single crystal permittivity tensor. The second-order forward bounds cannot be compared to this sea ice dataset, because they assume that the material is geometrically isotropic. This is not the case for columnar sea ice. However, the second-order bounds are applicable to granular ice.

As displayed in figure 2*a*, the polycrystalline bounds provide a much tighter bound than those bounding the permittivity of a general two-component material and statistically isotropic two-component material for sea ice. This makes sense because we are essentially applying a two-scale homogenization and including the additional information about rotation statistics. If we zoom in on the new polycrystalline bounds, we see that except for one data point, the bounds accurately capture the data for the corresponding volume fraction *ϕ* for the ice column (figure 3*a*). We suspect that the data point outside the bounds is explained by the variations in single crystal permittivities. However, we suggest that the three tighter bounds in figure 2*a* and the bounds in figure 3*a* might be viewed as bounds for a set of columnar sea ice permittivities for the whole column, where brine volume fractions *ϕ* vary in the interval 3.3%≤*ϕ*≤4.1%.

Although the bounds accurately capture the data, due the large variability and potential noise in both the data in [52] and the single crystal permittivity tensor obtained by numerical simulations, we further average the data in [52] and compare it to bounds corresponding to the single crystal permittivity tensor of the averaged value of *ϕ*=3.75%. This is displayed in figures 2*b* and 3*b*, which show that the bounds accurately capture the data point.

We further apply inverse polycrystalline bounds method to evaluate crystal orientation in sea ice samples using measured data of the ice permittivity. The goal of the inverse bounds is to estimate the mean crystal orientation potentially revealing the type of ice. For the data falling within a certain range of the brine fraction *ϕ*, the method gives bounds for the mean crystal orientation. The results are displayed in figure 4. The mean crystal orientation of the ice crystals in the ice column from the Arcone *et al.* experiments is between 8° and 30° off the vertical axis. Therefore, the ice is certainly columnar. We also examined images [52] of the typical sea ice structure, representative of the measured ice samples, and estimated that the mean crystal orientation should be between 11.5° and 19° off the vertical axis. These estimates are within the range of the inverse bounds on the mean single crystal orientation displayed in figure 4.

To compare the inverse mean crystal orientation bounds for columnar ice with the bounds for isotropic granular ice, we use an analytic model of a two-dimensional statistically isotropic polycrystalline material to calculate the effective permittivity of an isotropic polycrystal. From [23], it follows that the effective permittivity of such a polycrystalline composite coincides with the effective permittivity *θ*_{l} and the upper inverse bounds *θ*_{u} for the angle of deviation from the vertical axis are shown for 10 different effective complex permittivity values. The calculated bounds are in excellent agreement with the exact value of the deviation angle equal to 45°. The results demonstrate a significant difference in the reconstructed bounds for the mean orientation of a single crystal in columnar and in granular ice, which provides a foundation for distinguishing the types of ice using electromagnetic measurements.

## 5. Conclusion

We have developed both first- and second-order forward bounds on the effective complex permittivity ** ε*** for a polycrystalline material using the analytic continuation method. Additionally, we have derived first-order inverse bounds on the mean single crystal orientation for a polycrystalline material. The first-order forward bounds assume

*a priori*knowledge about the complex permittivity tensor for a single crystal and the mean single crystal orientation to bound

***. The second-order polycrystalline forward bounds further require the material to be geometrically isotropic in the polycrystalline Hashin–Shtrikman sense. The inverse bounds for the polycrystalline material assume knowledge of the effective permittivity**

*ε**** and the complex permittivity tensor for a single crystal to provide bounds for the mean crystal orientation. Comparison of the derived bounds with actual sea ice data show excellent agreement. These results provide a foundation for determining ice type with remote sensing techniques.**

*ε*## Funding statement

We gratefully acknowledge support from the Division of Mathematical Sciences and the Division of Polar Programs at the U.S. National Science Foundation (NSF) through grant nos. DMS-1009704, ARC-0934721, DMS-0940249, DMS-0602219 and DMS-1413454. We are also grateful for support from the Office of Naval Research (ONR) through grant no. N00014-13-10291. Finally, we would like to thank the NSF Math Climate Research Network (MCRN) for their support of this work.

- Received September 16, 2014.
- Accepted December 1, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.