Analytic Solutions to Reflection-Transmission Problem of Interface in Anisotropic Ice Sheet

+e 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. +en, 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. +is 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.


Introduction
e polar ice sheets play an important role in the global climate system. eir evolution and stability have great impact on climate change and sea level rising [1,2]. e 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]. e preferred fabrics and textures of ice crystals significantly affect the flow and evolution of ice [6][7][8].
erefore, 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 [9][10][11][12].
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]. ese 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 [22][23][24][25][26]. Additionally, some numerical methods have been developed for modeling the propagation of electromagnetic waves in the ice sheet. However, only a few [27][28][29] incorporated the anisotropy of the ice. Polarimetric radar is considered as the most potential tool for revealing the evolution of the anisotropic ice sheet. e major obstacle of using it in polar research is the limited understanding of radar wave propagation in the anisotropic ice sheet rather than hardware. is limited understanding deeply bottlenecks the data processing and interpretation of the polarimetric radar. e ice sheets can be simplified as a horizontally homogeneous stratified model with multiple anisotropic layers (isochronous layers) [27][28][29], which possess variable fabrics and permittivity tensors changed with the paleoclimate. erefore, 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 [30][31][32][33][34][35]. 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. e 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.
ere are some approaches to determine the reflections and transmissions of the EM waves at the interface of anisotropic media [39][40][41][42]. 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.

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. e 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. e 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.
e 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. e 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. e 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 secondorder permittivity tensor, which has a diagonal form along the principal axis direction: e fabric can be categorized based on the signatures of differentiation of permittivity along principal axes. e 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]. e four main fabric types can be categorized through the permittivity tensor: a random ice (Figure 1(a)) has permittivity elements of ε 1 � ε 2 � ε 3 ; a perfect single-pole ice (Figure 1(b)) has permittivity elements of ε 1 � ε 2 � ε ⊥ and ε 3 � ε // ; a vertical girdle fabric (Figure 1(d)) can be achieved when the permittivity elements satisfy ε 1 � ε // and ε 2 � ε 3 � (ε // + ε ⊥ )/2; and the permittivity tensor of an elongated single pole (Figure 1(c)) has ε 1 < ε 2 ≪ ε 3 .
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 x ′ denotes the semimajor axis of the indicatrix ellipses and y ′ 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 φ around z/z' axis.
Applying the coordinate transformation to the diagonal permittivity tensor (formula (1)) yields the generic form of the permittivity tensor in the measurement coordinate: 2 International Journal of Antennas and Propagation

Eigenanalysis of EM Wave Propagations in Anisotropic
Ice. e EM wave propagation in anisotropic ice-sheet layers can be described by differential-form Maxwell's equations in time-harmonic field: where E ⇀ i and H ⇀ i are the electric and magnetic fields in the ith ice sheet layer, ω is the angular frequency, μ 0 is the isotropic permeability of the i-th layer, and ε i 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 vector yields the matrix form of Maxwell's equations [43][44][45]: where e entries of Γ i 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 letting which yields an eigenvalue problem with the eigenvalues λ: In general, linear system (8) has four eigenvalues. erefore, the general solution of (7) can be written as the linear combination of eigenvectors of (8): or equivalently in matrix-vector notation: where a � a 1 a 2 a 3 a 4 , We can further reformulate (8) into where Equation (12) has two fundamental solutions: International Journal of Antennas and Propagation 3 For the normal incident EM wave along the -z/z′ axis, the wave field can be decomposed into waves I and II, where waves I and II are in the direction of x ′ and y ′ axis, respectively. In the measurement coordinate, waves I and wave II have rotation angles of φ and φ + π/2, respectively. e eigenvectors of the wave field in the i-th layer have two forms: International Journal of Antennas and Propagation where are intrinsic impedances of waves I and II in the i-th layer, subscripts 1 and 2 represent waves I and II, respectively, and φ i 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 e inverse matrices (18) and (19) resemble the matrices (15) and (16), except for a transpose operator, inverse of wave impedance, and a scaling factor. erefore, the inverse matrices (18) and (19) can be easily written once the dielectric tensor and azimuth angle are given. e 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 matrix and its inverse matrix (20) can be simplified into a matrix product of (15) and (22): International Journal of Antennas and Propagation e 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):

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. e 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. e reflection and transmission at the interface of the anisotropic half-space satisfy the following relation: We do matrix operation on equation (27) and get where R � R 11 R 12 Formula (28) illustrates that the eigenvector matrix of the first layer Na 1 and the inverse eigenvector matrix of the second layer Na Since the eigenvector matrix in each layer has two representations, there are four combinations to construct the matrix A � Na − 1 2 · Na 1 . We first consider a form with (20) and (25): Substituting (24) and (25) into (30) yields where Decomposing the eigenvector matrix into a diagonal matrix and an elementary matrix yields where I is the identity matrix, and b erefore, where and One can simplify (31) by substituting (36), (38), and (39): where Substituting (42) into (28) yields where one can solve for the reflection and transmission coefficients R and T: However, the expressions of R and T involve inverse matrix A − 1 a , which is typically difficult to obtain. Decomposing the matrix A 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 B a can be obtained from (36): Let a � (η 21 /η 11 ), b � (η 21 /η 12 ), c � (η 22 /η 11 ), d � (η 22 /η 12 ), and International Journal of Antennas and Propagation We can simplify the reflection coefficient (47) using (51): Similarly, substituting into (48) simplifies the transmission coefficient: For the other three combinations of A � Na − 1 2 · Na 1 , we can obtain the similar expression for the reflection and transmission coefficients using ctg(φ 2 − φ 1 ). For the special case of Δφ � φ 2 − φ 1 approximating kπ + π/2, the value of tg(φ 2 − φ 1 ) tends to infinity. So, we need to select appropriate combination for avoiding the infinity. where

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 International Journal of Antennas and Propagation Notice that the eigenvectors in the isotropic ice sheet only depend on the deflection angle at 0 and π/2. erefore,    Figure 7: e electric field components of the incident wave and reflection wave in the isotropic and anisotropic interface. Blue ellipse represents the indicatrix ellipse, and its long axis points to the direction of I; wave and short axis points to ; wave. Just incident wave of E H is chosen to demonstrate here, incident wave E V in axis y can be decomposed with the same method. the reflection and transmission coefficients (52) and (54) can be simplified to

Numerical Results
It had been mentioned that echo difference exists in radar exploration while rotating the antenna horizontally [19,[22][23][24]. erefore, we consider two numerical models to verify our aforementioned formulae for computing the reflection and transmission coefficients. e 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 ε ⊥ � 3.152 and ε // � 3.189 for model parameters, which had precisely been measured in the laboratory [20,21].

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 ε 1 � 3.152, ε 2 � 3.189, and ε 3 � 3.189. In model 1, we fix media coordinates for both layers and keep Δφ � φ 2 − φ 1 � π/6, where φ 1 and φ 2 are azimuth angles of the first and second layer relative to the measurement coordinate. en, we rotate the measurement coordinate from 0 to 2π and calculate each reflection and transmission coefficient at each rotating angle. is operation emulates the polarizing antenna rotates around the horizontal plane inversely. e 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 2π 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 Δφ � φ 2 − φ 1 from 0 to 2π, where φ 1 and φ 2 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.

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 ε 1 � ε 2 � ε 3 � 1 as the free air layer, and the lower layer is anisotropic with ε 1 � 3.152, ε 2 � 3.189, and ε 3 � 3.189 as the real ice layer. In this model, we fix media coordinates for both layers and keep φ 2 − φ 1 � 0 and φ 2 − φ 1 � − π/6 ( Figure 5), where φ 1 and φ 2 are azimuth angles of ice fabric in the first and second layer. We rotate the measurement coordinate from 0 to 2π 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 − π/6, while the difference of azimuth angle changes from 0 to − π/6. is can be used for tracking the azimuth of the maximum echo in radar survey on the ice sheet.

Physical Meaning of Eigenvalues and Eigenvectors.
e 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: e column vectors of eigenvector matrices in (15) Equation (63) shows that the ratios E x /E y and H x /H y are only functions of the media angle φ, which is measured between the measurement and the media coordinate. e 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 φ. erefore, 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 φ.

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 timeshare 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 φ 2 relative to the measurement coordinate. E H is the incident wave in the isotropic layer. It can be decomposed into E HI and E HII along the semimajor and semiminor axes of the indicatrix ellipse. ese two incident waves generate two reflection waves E R HI and E R HII from the interface:

Data Availability
e data used to support the findings of this study are included within the article.

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. e angular separation in Model 31 is zero. e angular separation in Model 32 is set -π/6, which makes its curves shift right π/6. (Supplementary Materials)