Abstract

The rheology and evolution of the polar ice sheet are deeply influenced by the anisotropy of ice crystals. Studying the anisotropy of ice crystals can help to well understand and predict the behavior of the polar ice sheet and then the sea level rising and global climate change. In this paper, firstly, we deduce the expression of eigenvalues and eigenvectors of anisotropic media, which are determined by permittivity tensor and geometry of media. Then, the analytic formulas of reflection and transmission coefficients are derived directly by matrix transformation. Some models with real ice parameters are tested, and they present some special features at the anisotropic interface. We also discuss the physical meanings of eigenvalues and eigenvectors and the geometry analyzing to polarimetric radar. This analytic solution reveals the functional relationship between the macroradar reflection and the microphysical properties of ice crystals, which provides a feasibility of ice fabric identification by polarimetric radar detection.

1. Introduction

The polar ice sheets play an important role in the global climate system. Their evolution and stability have great impact on climate change and sea level rising [1, 2]. The huge thick ice sheet [3] (more than 2,000 m on average) is composed of numerous microice crystals. A single ice crystal shows noticeable anisotropy in its response to stress, specified by its exclusive C-axis, which is orthogonal to its basal plane. However, the bulk of ice crystals may exhibit preferred crystal orientations (fabric) and variable shapes and sizes (texture) under the giant stress in the ice sheet [4, 5]. The preferred fabrics and textures of ice crystals significantly affect the flow and evolution of ice [68]. Therefore, the fabric alignment in the ice sheet is a key indicator for understanding the behavior of the ice sheet. Unfortunately, this information can only be obtained from very few deep ice cores in Antarctic and Greenland [912].

Glen [13] presented an isotropic ice flow law, which successfully explained the early field observations and had been widely accepted. However, recent field measurements showed some discrepancies with the results predicted by Glen’s flow law, given the improvements from equipment detection precision and extent [14]. These discrepancies implied that the anisotropy of the ice crystal could not be ignored in understanding the ice sheet evolution. Considering the importance of ice anisotropy in the ice sheet evolution [9, 14], how to recognize the fabric and its distribution in the ice sheet and how to understand its anisotropic behavior are still the leading and challenging fields in polar research [15].

RES (or radio-echo sounding) has been widely employed as the standard facility in ice sheet expeditions since the 1950s. It achieves remarkable success in the polar exploration due to its high efficiency and accuracy [16, 17]. Early works since the 1950s had reported the elliptical polarization in RES echoes in ice sheet survey [18]. Hargreaves [19] suggested birefringence as the cause of echo difference along different polarized orientations and proposed an anisotropic expression of the permittivity with tensor. It is till the late 1990s that the permittivity tensor of the ice crystal within the frequency band of radar waves is precisely measured in the laboratory [20, 21]. Some field surveys have been conducted to address the relation between the polarization of the EM wave and the ice fabric in the ice sheet [2226]. Additionally, some numerical methods have been developed for modeling the propagation of electromagnetic waves in the ice sheet. However, only a few [2729] incorporated the anisotropy of the ice. Polarimetric radar is considered as the most potential tool for revealing the evolution of the anisotropic ice sheet. The major obstacle of using it in polar research is the limited understanding of radar wave propagation in the anisotropic ice sheet rather than hardware. This limited understanding deeply bottlenecks the data processing and interpretation of the polarimetric radar.

The ice sheets can be simplified as a horizontally homogeneous stratified model with multiple anisotropic layers (isochronous layers) [2729], which possess variable fabrics and permittivity tensors changed with the paleoclimate. Therefore, while considering the propagation of the EM wave in the ice sheet, the solution for this problem can be categorized into two sequential steps: the reflection and transmission at the anisotropic interface and the propagation in stratified and anisotropic media. For stratified anisotropic media, Berreman’s 4 × 4 propagation matrix is the most popular method used in optic and other EM fields [3035]. Ursin [36] presented a propagation matrix method commonly used for elastic and EM wave propagating in horizontally layered media. Sluijter et al. [37] used the polarized ray-tracing method for analyzing the inhomogeneous anisotropic media. However, in ground penetrating radar, the wave sources are often transient pulses, not the same as the continuous wave used in EM and optic fields. The echoes of radar can be considered as the convolution response of a ray series (or geometrical optics series) [38] with the excitation wavelet. Each ray series is the combination of multiple reflections and transmissions at interfaces along the ray path. So, the reflection and transmission at the interface is a key basis for understanding the EM propagation in stratified media. There are some approaches to determine the reflections and transmissions of the EM waves at the interface of anisotropic media [3942]. In this paper, we present a matrix transformation method to derive the direct expressions of reflection and transmission at the interface. It will be useful for the propagation in stratified and anisotropic media.

2. Methodology

2.1. Simplification of the Electromagnetic Model

RES survey and deep ice core analysis have proved that, in small scale, the deep ice sheet can be simplified as horizontal stratified ice layers. The ice sheet is often composed of four main fabric types [6, 27].

Figure 1 demonstrates the Schmidt diagrams with the four fabrics in the ice sheet. C-axis is a direction orthogonal to the crystal basal plane and often acts as the unique indicator of a crystal. The gray points or filled areas in Figure 1 indicate the possible C-axis directions or the gathering shape of C-axis (fabric) in polycrystalline ice. The symbols x and y indicate the principal axes in the horizontal plane, and the z-axis is perpendicular to the horizontal plane. Random ice (Figure 1(a)) often appears in the shallow part of the ice sheet because the stress is faint, and the C-axis points any direction kept while snow fell on the surface. The perfect single pole (Figure 1(b)) and vertical girdle (Figure 1(d)) are the two extreme states of fabric with ice crystal deformation under the stress and shear. The elongated single pole (Figure 1(c)) is the transient state between the perfect single pole and the vertical girdle under shear. Mathematically, anisotropic crystals in the ice sheet can be described by the second-order permittivity tensor, which has a diagonal form along the principal axis direction:

The fabric can be categorized based on the signatures of differentiation of permittivity along principal axes. The permittivity of a crystal along or perpendicular to its C-axis is denoted as and , respectively, which had precisely been measured in the laboratory [20, 21]. The four main fabric types can be categorized through the permittivity tensor: a random ice (Figure 1(a)) has permittivity elements of ; a perfect single-pole ice (Figure 1(b)) has permittivity elements of and ; a vertical girdle fabric (Figure 1(d)) can be achieved when the permittivity elements satisfy and ; and the permittivity tensor of an elongated single pole (Figure 1(c)) has .

In radar analysis, EM wave propagations are typically described in two coordinates: measurement and media coordinates. Measurement coordinate is used to describe the amplitude and propagation direction of EM components. Media coordinate is used to define the principal axes of permittivity tensor and geometry of media. Figure 2 presents the relative position of the measurement and media coordinates, where denotes the semimajor axis of the indicatrix ellipses and denotes the semiminor axis. Note that the z-axis in the measurement coordinate coincides with the z’-axis in the media coordinate because EM waves propagate as plane waves with normal incident angle into the ice sheet, where isochrones are generally horizontal. In the layered ice model, the measurement coordinate is shared among all layers, whereas the media coordinate differs in each layer and is simply a rotation of the measurement coordinate with an azimuth angle .

Applying the coordinate transformation to the diagonal permittivity tensor (formula (1)) yields the generic form of the permittivity tensor in the measurement coordinate:

2.2. Eigenanalysis of EM Wave Propagations in Anisotropic Ice

The EM wave propagation in anisotropic ice-sheet layers can be described by differential-form Maxwell’s equations in time-harmonic field:where and are the electric and magnetic fields in the i-th ice sheet layer, is the angular frequency, is the isotropic permeability of the i-th layer, and is the permittivity tensor of the i-th layer. In the measurement coordinate, we suppose z-axis is the direction of EM propagation. It satisfies the case of normal incidence. So, electric field E and magnetic field H just exist in the plane of x- and y-axis determining. Combining electric and magnetic fields into one vectoryields the matrix form of Maxwell’s equations [4345]:where

The entries of are expressed in terms of the azimuth angle and permittivity tensor in the i-th layer. For vertically propagating plane waves with wave number , the eigensolutions of (5) can be found by lettingwhich yields an eigenvalue problem with the eigenvalues :

In general, linear system (8) has four eigenvalues. Therefore, the general solution of (7) can be written as the linear combination of eigenvectors of (8):or equivalently in matrix-vector notation:where

We can further reformulate (8) intowhere

Equation (12) has two fundamental solutions:

For the normal incident EM wave along the - axis, the wave field can be decomposed into waves I and II, where waves I and II are in the direction of and axis, respectively. In the measurement coordinate, waves I and wave II have rotation angles of and , respectively. The eigenvectors of the wave field in the i-th layer have two forms:whereare intrinsic impedances of waves I and II in the i-th layer, subscripts 1 and 2 represent waves I and II, respectively, and is the azimuthal rotation angle of the i-th layer in the measurement coordinate with respect to the media coordinate. Subscripts (1) and (2) indicate the two forms of matrices in (15) and (16), which are eigenvectors composed of tangent and cotangent functions of the azimuth angles, respectively. Taking the inverse of matrices (15) and (16) yields

The inverse matrices (18) and (19) resemble the matrices (15) and (16), except for a transpose operator, inverse of wave impedance, and a scaling factor. Therefore, the inverse matrices (18) and (19) can be easily written once the dielectric tensor and azimuth angle are given.

The calculation of reflection-transmission coefficients needs normalized eigenvectors, so we first normalize each column of the matrix in (15) by the L2 norm of the column vector:where

We introduce a diagonal matrixand its inverse matrix

(20) can be simplified into a matrix product of (15) and (22):

The corresponding inverse matrix of (20) can be simplified using (18) and (23):

We can also write (24) and (25) in terms of the cotangent functions of rotation angles from (16) and (19):

2.3. Reflection and Transmission Coefficients for the Anisotropic Interface in the Ice Sheet

We use the 2 × 2 reflection matrix R and the 2 × 2 transmission matrix T for the problem of EM propagation at the anisotropic interface. Considering an interface between two anisotropic ice-sheet layers, the reflection matrix R can be defined to relate the reflected upgoing waves to the downgoing incident waves. The transmission matrix T relates the transmitted downgoing waves to the downgoing incident waves. Above the interface, the total wave field consists of the incident downgoing wave and the reflected upgoing wave. Below the interface, the total wave field has only downgoing transmitted waves. The reflection and transmission at the interface of the anisotropic half-space satisfy the following relation:

We do matrix operation on equation (27) and getwhere

Formula (28) illustrates that the eigenvector matrix of the first layer and the inverse eigenvector matrix of the second layer are necessary for calculating and . Since the eigenvector matrix in each layer has two representations, there are four combinations to construct the matrix . We first consider a form with (20) and (25):

Substituting (24) and (25) into (30) yieldswhere

Decomposing the eigenvector matrix into a diagonal matrix and an elementary matrix yieldswhere

I is the identity matrix, and and are block diagonal matrices. Since the product of two block diagonal matrices is still block diagonal, we arrive at

Therefore,where

Let

and

One can simplify (31) by substituting (36), (38), and (39):where

Substituting (42) into (28) yieldswhere one can solve for the reflection and transmission coefficients R and T:

However, the expressions of R and T involve inverse matrix , which is typically difficult to obtain. Decomposing the matrix into block diagonal matrices and elementary matrices as shown in (42) and substituting them into (45) and (46) yield

In this way, complex calculation of R and T can be decomposed into product or sum between the matrix on the diagonal matrix F and the submatrix of B and its inverse matrix, where the inverse of matrix can be obtained from (36):

Let , and

We can simplify the reflection coefficient (47) using (51):

Similarly, substitutinginto (48) simplifies the transmission coefficient:where

For the other three combinations of , we can obtain the similar expression for the reflection and transmission coefficients using . For the special case of approximating , the value of tends to infinity. So, we need to select appropriate combination for avoiding the infinity.where

2.4. Special Case: Isotropic and Anisotropic Interface

We now consider the special case of isotropic ice media, where the permittivity tensor is isotropic, and (6) can be simplified to

Following the aforementioned derivations, the corresponding eigenvalues and eigenvectors are

Notice that the eigenvectors in the isotropic ice sheet only depend on the deflection angle at 0 and . Therefore, the reflection and transmission coefficients (52) and (54) can be simplified to

3. Numerical Results

It had been mentioned that echo difference exists in radar exploration while rotating the antenna horizontally [19, 2224]. Therefore, we consider two numerical models to verify our aforementioned formulae for computing the reflection and transmission coefficients. The aforementioned formulae have no special limit to the range of permittivity. However, for test and verification of the real reflections in the ice sheet, we use the permittivities and for model parameters, which had precisely been measured in the laboratory [20, 21].

3.1. Reflection and Transmission Coefficients for an Interface between Anisotropic Layers

Firstly, we construct a two-layer ice model, where both layers are anisotropic and the values of the entry permittivity tensors are , , and . In model 1, we fix media coordinates for both layers and keep , where and are azimuth angles of the first and second layer relative to the measurement coordinate. Then, we rotate the measurement coordinate from 0 to and calculate each reflection and transmission coefficient at each rotating angle. This operation emulates the polarizing antenna rotates around the horizontal plane inversely.

The reflection and transmission coefficients are drawn in Figure 3. Both reflection and transmission coefficients are kept constant, while rotating the measurement coordinate from 0 to in the horizontal plane. When we change the permittivity or geometry , the values of the reflection and transmission coefficients will change too, but they still remain constant in a new level while measurement coordinate rotating.

In model 2, we fix the measurement coordinate and scan the difference of media azimuth angle from 0 to , where and are azimuth angles of the first and second layer. In this case, we observe that the reflection and transmission coefficients vary periodically in Figure 4.

3.2. Reflection and Transmission Coefficients for an Interface between Isotropic and Anisotropic Layers

Considering the real survey environment of the ice sheet, we construct another two-layer ice model, in which the upper layer is isotropic with as the free air layer, and the lower layer is anisotropic with , , and as the real ice layer. In this model, we fix media coordinates for both layers and keep and (Figure 5), where and are azimuth angles of ice fabric in the first and second layer. We rotate the measurement coordinate from 0 to in the horizontal plane.

In this case, we observe that reflection and transmission coefficients vary periodically in Figure 5 while we rotate the measurement coordinates. Comparing the variation of reflection and transmission, we can see the curves shift right , while the difference of azimuth angle changes from 0 to . This can be used for tracking the azimuth of the maximum echo in radar survey on the ice sheet.

4. Conclusion and Discussion

4.1. Physical Meaning of Eigenvalues and Eigenvectors

The eigenvalues and eigenvectors contain information about the relation between electromagnetic components and physical properties of media when the EM wave propagates in anisotropic media. Multiplying (13) by the angular frequency yields the upgoing and downgoing wave vectors of I and II waves:

The column vectors of eigenvector matrices in (15) and (16) show the proportional relation of EM components, . The four columns in eigenvector matrices represent the normalized component of the upgoing and downgoing wave vectors corresponding to I and II waves, respectively. In anisotropic media,

Equation (63) shows that the ratios and are only functions of the media angle , which is measured between the measurement and the media coordinate. The intrinsic impedance of media relates the magnetic and electric fields, as shown in (64). So, the eigenvector matrix describes the proportion relationship between EM components of I and II waves and proves that eigenvectors are only determined by intrinsic impedance of media and the deflection angle . Therefore, we can conclude that the eigenvector matrix describes the proportional relation of EM components of I and II waves. Moreover, the eigenvectors are only determined by the intrinsic impedance of media and the deflection angle .

4.2. Simplification of the Formulae for Reflection and Transmission Coefficients

For natural media, and in formulas (40) and (41) tend to the unit matrix. In formula (52),

In formula (54), F1, F2, F3, and F4 also tend to 1, so formulas (52) and (54) can be simplified to

4.3. Polarimetric Radar Survey and Geometry Analysis

A typical configuration of polarimetric radar consists of two orthogonal transmitters and two orthogonal receivers. Four transmitter-receiver pairs are possible, either copolarized or cross-polarized, which are shown in Figure 6. Conventional polarimetric radar implements configuration with time-share single transmitter to double receiver alternately. Orthogonal echo differences can be used to infer the anisotropic feature of the media.

For an interface between isotropic-anisotropic layers, we can assume that the measurement coordinate is {X, Y} in Figure 7, and the media coordinate of the anisotropic layer has deflection angle relative to the measurement coordinate. EH is the incident wave in the isotropic layer. It can be decomposed into EHI and EHII along the semimajor and semiminor axes of the indicatrix ellipse. These two incident waves generate two reflection waves ERHI and ERHII from the interface:

Therefore, the reflection from the interface along the x-axis is the vector sum of ERHI and ERHII along the x-axis. The reflection coefficient R11 can be calculated asand similarly, the reflection coefficient R12 can be computed as

For incident wave along the y-axis, we can derive the reflection coefficients R21 and R22 in a similar fashion, which have been shown in (60).

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was based on the polar expeditions of China. The authors thank all members in CHINARE 32, especially Dr. Guo Jinxue and Dr. Cui Xiangbing. They also thank Prof. Tian Gang in the School of Earth Science of Zhejiang University for his useful suggestion. This work was supported by the National Natural Science Foundation of China under grants 40874060, 41176167, and 41676181, the National Basic Research Program of China (973 Program) under grant 2012CB957702, and the Fundamental Research Funds for the Central Universities under grant 2019XZZX005-2-03.

Supplementary Materials

In the directory of data, there are four excel files in it. Every excel file consists of two parts: the first part is the model parameter setting of the layers and includes 3 relative permittivity and azimuth angle Fi in this layer, and the second part is the calculated result according to the model and includes 4 reflection coefficients and 4 transmission coefficients calculated while the antenna rotated around 360 degrees in the horizontal plane. Model0.xlsx: a model and its calculated results for comparison with the past traditional isotropic formula of reflection and transmission coefficients. For the convenience of rapid verification, the relative permittivity is set with simple 1 and 4. No figure for this model is used in paper. Model1 for Figure 3.xlsx: a model with two anisotropic layers and calculating results for Figure 3. In this model, the angular separation between the up and down layers is fixed and the direction of antenna rotate is around 360 degrees in the horizontal plane. Model2 for Figure 4.xlsx: a model with two anisotropic layers same as model 1 and the calculated results for Figure 4. In this model, the azimuth angle of the up layer is fixed, and the azimuth angle of the down layer changes with the antenna rotating around 360 degrees in the horizontal plane. So, the angular separation keeps changing with antenna rotation. Model31 for Figure 5.xlsx and Model32 for Figure 5.xlsx: two models with two different angular separation layers and calculated results for Figure 5. In these models, the up layer is an isotropic layer and the bottom layer is an anisotropic layer. The angular separation in Model 31 is zero. The angular separation in Model 32 is set –π/6, which makes its curves shift right π/6. (Supplementary Materials)