Remote Sensing of the Polar Ice Zones with HF Radar

Abstract

Radars operating in the HF band are widely used for over-the-horizon remote sensing of ocean surface conditions, ionospheric studies and the monitoring of ship and aircraft traffic. Several hundreds of such radars are in operation, yet only a handful of experiments have been conducted to assess the prospect of utilizing this technology for the remote sensing of sea ice. Even then, the measurements carried out have addressed only the most basic questions: is there ice present, and can we measure its drift? Recently the theory that describes HF scattering from the dynamic sea surface was extended to handle situations where an ice cover is present. With this new tool, it becomes feasible to interpret the corresponding radar echoes in terms of the structural, mechanical, and electrical properties of the ice field. In this paper we look briefly at ice sensing from space-borne sensors before showing how the persistent and synoptic wide area surveillance capabilities of HF radar offer an alternative. The dispersion relations of different forms of sea ice are examined and used in a modified implementation of the electromagnetic scattering theory employed in HF radar oceanography to compute the corresponding radar signatures. Previous and present-day HF radar deployments at high latitudes are reviewed, noting the physical and technical challenges that confront the implementation of an operational HF radar in its ice monitoring capability.

1. Introduction

The effects of climate change are nowhere more evident than in the Arctic, where the inter-annual retreat of the ice cover is occurring at a rate exceeding that predicted by almost all major climate models [1,2]. Moreover, the reduction in ice extent is much more pronounced around the extent minimum, normally in September, than it is near the coverage maximum in February–March. Figure 1 illustrates this seasonal distinction with the latest data from the two periods. In this figure, adapted from [3], the magenta lines mark the median boundaries over the period 1981–2010.

Apart from its impact on global warming via several feedback mechanisms, the reduced ice cover is opening up a range of human activities hitherto either impossible or prohibitively too expensive to consider seriously. Foremost among these is the prospect of a massive expansion of shipping through the Northern Sea Route, taking advantage of the change in navigability, primarily during the summer months, when the ice conditions are highly dynamic and wind–wave activity in the Arctic increases, with greater fetch available for wave growth. Other activities that are likely to increase include extraction of undersea oil, gas and mineral resources, expanded fishing and harvesting of other aquatic biota, and establishment of permanent facilities typically linked with claims of national sovereignty. In all these domains, and for the over-arching mission of climate monitoring and modelling, a key requirement is detailed, synoptic, and timely information about the sea ice; and not merely its distribution but its physical characteristics.

Ideally, one might wish to deploy a dense network of in situ sensors that are able to measure all the parameters of concern, but the impracticality of this option obliges us to a system consisting of fewer sensors. In order to achieve wide area surveillance, an individual sensor must have to be able to interrogate an entire region of interest from one fixed location, or it must travel across the region, sampling as it goes. The latter case is the prevailing technology–space-borne sensors programmed with orbital parameters constrained to meet the ground coverage needs. Despite the impressive capabilities of space-based sensors, the technologies presently employed to monitor sea ice fall short of what is needed, so in this paper we shall explore the prospect of a fixed site sensor with true wide area surveillance capability, namely HF ‘over-the-horizon’ radar.

In the following section, we review the ice monitoring technologies deployed in space-based systems, and note their limitations as well as their capabilities. Then, in Section 3, we turn our attention to HF radar, emphasizing its reliance on the geophysical environment to support the special propagation modes that enable it to observe the Earth’s surface out to ranges of several thousands of kilometers via skywave and several hundreds of kilometers via surface wave. The scattering theory used to model the reflection of HF radio waves from the sea surface is reviewed and generalized to accommodate surfaces with different dynamical behavior, as is the case when ice is present. Armed with this understanding, we review previous experiments that use HF radar at high latitudes, where echoes from sea ice have been observed, and identify the attribute that opens the way to the prospect of sea ice characterization.

Having introduced our sensor, in Section 4 we proceed to survey the properties of the target of interest: sea ice in its various forms. The evolution of ice fields is governed by a host of meteorological and oceanographic phenomena, but the accumulated studies of many researchers lead us to the conclusion that models based on a reduced set of parameters can describe the rheology of ice in nearly all these forms. The key feature of interest from the HF radar perspective is the dispersion relation of those waves in the ice field which have been excited by surface gravity waves from the open sea penetrating the ice field. It is tempting to think of these excited waves as merely the original open sea waves subject to attenuation due to viscous and frictional losses, but this interpretation cannot be pushed far because a whole menagerie of oscillatory behavior can appear. To elucidate the ways in which models attempt to represent the dynamical properties of sea ice, we derive the dispersion relations for simple models and sketch the directions taken by the more sophisticated approaches, which add increasing complexity to the ice mechanics. Studies of model validation are reviewed as a guide to model selection and assignment.

Next, in Section 5, we describe the application of the scattering theory of Section 3 to the ice models of Section 4. The computational framework is presented, and examples of HF signatures of sea ice are shown and used to examine the questions of ice species identification and ice parameter measurement accuracy. The discussion to this point focusses on local properties and signatures, but in Section 6 we extend the perspective to consider large-scale attributes of ice fields and the application of the action balance equation to model the spatial evolution of the ice field and hence its radar signatures. The prospect of a reciprocal benefit, wherein HF radar could enhance our understanding of the source terms in the action balance equation, is identified.

Section 7 looks at the practical issues involved with HF radar operations in the ice zone, especially factors that limit radar capability, together with a discussion of measures that can mitigate those factors. Section 8 surveys the HF radars presently in operation or under development that could be employed to explore the ideas proposed in this paper.

In the Conclusion, we detail our assessments of the capabilities that HF radar signatures offer.

2. Space-Borne Ice Sensors and Their Observables

The greatest spatial coverage is afforded by space-based sensors, which we can separate into passive and active classes.

Passive microwave radiometers (PMRs) have been a mainstay for many years [4]. Apart from the relatively coarse spatial resolution, with typical cell dimension ~25 km, passive microwave retrievals experience difficulty when the sea ice is undergoing surface or lateral melt processes, as well as being degraded somewhat by the presence of heavy cloud cover. In addition, the emissivity of sea ice is itself highly variable, apart from a systematic dependence on the type of ice. Thus, during summer, when we can expect greatest human activity to occur, PMRs are often compromised. Nevertheless, modern implementations, such as AMSR2 [4], employ multiple frequency bands and, when conditions are favorable, can deliver up to seven products, namely, sea ice concentration, integrated total water vapor, liquid water path, sea surface temperature (SST), ice surface temperature and multi-year ice fraction [5].

Optical and multispectral imagers such as OLCI (ocean and land color instrument) on Sentinel-3 [6] rely on the reflection of natural sources of illumination, and are therefore hindered by clouds and darkness. A key advantage of these instruments is that the images are readily interpretable to human users, who can draw on their experience with viewing from airborne platforms. In addition, with high spatial resolution ~10 m, they provide the ability to detect scene changes due to both ice field slow-time evolution and drift, but they are not able to directly measure dynamics.

Active space-borne sensors such as altimeters, scatterometers and synthetic aperture radar (SAR) exploit the sensitivity of the radar scattering mechanism to radar frequency, scattering geometry and polarization, as well as to ice characteristics such as surface geometry, thickness, topography and age (which influences the ice dielectric properties). In broad terms, these instruments fall into three categories.

Regarding altimeters, it has been asserted that ‘satellite altimetry is the only method to monitor global changes in sea-ice thickness and volume over decades’ [7]. Altimeters, such as that on CryoSat-2 [8], employ either microwave radars or lasers to measure the distance to the reflecting ice surface; combining this with an estimate of the local sea level yields an estimate of the ice freeboard, and hence thickness, but it can only be performed directly under the satellite’s path and is subject to bias from the presence of snow. As the measurements are only one-dimensional and have a revisit time of about 27 days, they have restricted applications as stand-alone observations but are invaluable as a validation or calibration system for other sensors. A recent observation relevant to the HF radar technique reported in this paper is the detection and interpretation of ice distortions due to penetrating long ocean waves [9].

Scatterometers (see [10] for a recent review) are conventional real aperture radars that measure the Earth surface echoes at some set of parameter values, especially polarization states, and use these to infer the nature of the reflecting surface and, in some cases, the underlying body. In the context of the cryosphere, scatterometers can detect the presence of ice and aspects of its dielectric properties, often enabling one to discriminate between first-year ice and multi-year ice. Swathe widths of ~500 km provide a reasonable coverage rate, at a spatial resolution of ~25 km, with revisits as frequent as daily, so crude estimates of ice drift can be derived.

Synthetic aperture imaging radar (see [11] for a recent review) exploits the high spatial resolution gained through aperture synthesis, typically 3–15 m, at the expense of various image distortion mechanisms and speckle. As with scatterometers, the exploitation of the polarization domain is central to the retrieval of detailed ice characteristics, but information on ice thickness must be inferred rather than measured directly [12]. In fact, there are several ways in which ice thickness can be estimated from SAR data, but the most effective techniques employ a static version of the dynamic concept that we propose to use at HF, namely, examination of the perturbations of the ice surface due to waves penetrating from the adjacent open sea. SAR achieves this in the spatial domain, by virtue of its high spatial resolution, revealing patterns in scattered intensity, which indicate the presence of long waves; this was first exploited in [13] to estimate ice properties by measuring the changes in wavelength as long waves entered an ice field, and testing this against model dispersion relations. Subsequently, more attention has been paid to changes in wave amplitude than wavelength, as this yields greater sensitivity. We will explain this later in the paper. These techniques have since been applied in numerous scenarios [14,15,16,17,18,19,20,21,22]. At present, the revisit period for individual operational SAR systems is quite long, typically 3–15 days, preventing meaningful estimation of drift and ice growth and decay. Constellations of satellites, such as the RADARSAT Constellation Mission (RCM), launched in June 2019, can achieve proportionately more frequent revisits. In the case of the RCM, three satellites are in orbit at an altitude of 600 km, with an orbit time of 96 min. Together, they are able to deliver daily revisits of 90% of the Earth’s surface.

Nominal capabilities of individual space-based sensors, as they relate to ice monitoring, are summarized in

Table 1. Observables and nominal measurement capabilities of space-borne ice sensors. Lower case letters indicate a partial capability.

Observable	Passive Sensors		Active Sensors
	Visible and IR Imagers	Microwave Radiometers	Altimeters	Scatter- Ometers	LEO Synthetic Aperture Radar
ice field edge	Y	Y	Y	Y	Y
ice type	Y	y	N	y	Y
ice thickness	N	Y	Y	N	y
floe size distribution	Y	N	y	N	Y
permittivity	y	Y	N	Y	y
ice drift	Y	Y	N	y	Y
ice dynamics	N	N	N	N	N
resolution (horiz.)	10–30 m	5–25 km	10–70 m	25 km	1–15 m
revisit time	1–2 days	1 day	27 days	7–14 days	3–15 days

.

3. HF Radar and Its Observables

3.1. Propagation Modes and System Characteristics

Unlike the microwaves employed by active space-borne sensors, with frequencies mainly in the range 1–35 GHz, radio waves in the HF band, 3–30 MHz, are able to propagate from a ground-based source to the terrestrial surface far beyond the visible horizon. There are several mechanisms that support this non-line-of-sight propagation; of these, two provide the modes exploited by HF radar:

• Skywave–reflection from the Earth’s ionosphere, at heights between 60 and 400 km;
• Surface wave–diffraction around the curved Earth’s surface. At long ranges this takes the form of a localized electromagnetic wave seemingly attached to the surface, accompanied by electrical currents driven in the conducting seawater.

The skywave mode provides illumination of the distant Earth surface to ranges of several thousand kilometers by a single reflection, and further when multiple bounces in the Earth-ionosphere waveguide are considered. The price one pays for this vast coverage is distortion of the radar signals due to the fluctuations in the ionospheric plasma: the ‘space weather’. The surface wave mode avoids the ionospheric distortion but the diffracted field strength decays much more rapidly, supporting radar operations out to ranges of several hundred kilometers. We illustrate the resulting coverage in Figure 2 for monostatic radars, that is, radars where the transmit and receive sites are co-located.

Unlike the space-borne sensors, constrained by their orbital parameters and correspondingly long revisit times, HF radars can provide continuous (or persistent) surveillance. In the case of HFSWR (HF surface wave radar), this is usually achieved by simultaneous analysis of the echoes from all the resolution cells within the illuminated footprint, perhaps several thousand, updated on a timescale of typically 50–600 s, depending on the type of observation. This time interval or dwell time determines the achievable resolution in the Doppler domain, and thus in the fineness of detail about the dynamics of the scatterers in the illuminated zone. In velocity terms, it is often possible to resolve scatterers whose speeds differ by less than 0.1 ms⁻¹. Skywave radars similarly analyze the echoes from the hundreds of resolution cells within the selected dwell illumination region (DIR), such as that shown in Figure 2, usually employing somewhat shorter dwell times, often in the range 1–200 s. If there is a need for monitoring a number of DIRs, each of an area of perhaps 200,000 km², the radar will scan them sequentially, in which case the revisit time may be extended to 10–15 min, which is still adequate to monitor even rapid changes in sea ice conditions.

To emphasize the importance of persistent observation, we note that ice floes in the Arctic can move surprisingly quickly, up to 40 km in a day [23], and the ice edge can advance or retreat over a wide front by hundreds of kilometers between revisits of a spaceborne sensor and can breakup under wave action at speeds up to 16 km/h.

As with any radar, the spatial resolution of HF radars is governed mainly by three factors—the radar signal frequency, the waveform bandwidth, and the receive antenna aperture employed—but the numerical values of these differ greatly from those of microwave radars. Resolution in range (time delay) is limited by access to clear channels in the HF spectrum and by propagation effects, but commonly lies between 1.5 and 6 km when remote sensing missions are undertaken. The cross-range dimension of a resolution cell is set by the aperture diffraction limit, so a 1 km antenna array, operating at a frequency of 10 MHz, yields nearly a 3 km cross-range resolution at 100 km range or 60 km at 2000 km. Thus, to cite broadly representative values, individual resolution cells may have areas in the range 4–20 km² in the case of HFSWR systems, while large skywave radars tend to operate with a resolution cell area of 25–100 km² when studying ocean surface conditions.

More details about HF radars can be found in [24,25,26].

It is evident from these numbers that HF radar has the attributes required to maintain a constant watch over wide areas of the terrestrial surface at a spatial resolution commensurate with the level of detail desired and with a Doppler resolution far superior to that achievable with space-borne microwave radars. The question at the core of our investigation is this: what can an HF radar tell us about sea ice? Researchers familiar with the difficulties of deploying sensors in the ice zone might also ask whether it is feasible to construct and operate HF radars at high latitudes where they can survey the regions of interest.

To begin to answer the first question, we must look at the well-established oceanographic applications of HF radar and, in particular, at the scattering mechanism that explains the varied radar signatures of ocean surface structure and dynamical processes.

3.2. Radio Wave Scattering from the Ocean Surface

It should be recalled that, except near vertical incidence, microwave radar echoes arise predominantly by Bragg scatter, so the structures on the ocean surface that contribute to the scattered field are the short capillary and capillary–gravity waves with wavelengths ranging from millimeters to a decimeter or two, depending on the radar frequency. These waves carry only a tiny part of the ocean surface energy; they ride on the much more energetic surface gravity waves and swell, with wavelengths of 5–200 m and amplitudes of the order of a meter. It follows that remote sensing at microwave frequencies entails the use of models that attempt to diagnose the distribution of the energetic gravity waves from the spatial texture of echoes generated by the ‘marker’ capillary–gravity waves.

The texture arises in part because the scattering process is geometry-dependent, so the varying slope of the underlying long waves impresses a modulation on the radar echoes. To complicate matters, the short marker waves are not uniformly distributed over the longer waves. In addition, the short waves are advected by the non-uniform surface flows associated with the long waves, and there is also a nonlinear hydrodynamic coupling mechanism that modulates the short waves. As a consequence, the simultaneous presence of highly complex electromagnetic and hydrodynamic processes forces us to retreat from a purely analytical approach and instead adopt hybrid models based on a mixture of theory and empirical methods.

In contrast, the wavelengths of radio waves in the HF band are comparable with those of the most energetic surface gravity waves, and much greater than the surface wave amplitudes except in very high sea states. Moreover, at HF, the intricacies of centimeter-scale capillary–gravity waves can be ignored; they contribute insignificantly to the overall scattered field. Thus, on both electromagnetic and hydrodynamic fronts, the physics that needs to be considered is far simpler than its microwave counterpart, so an analytical approach is feasible. The most widely used approach, the small perturbation method (SPM) [27], solves the full radar scattering problem as a power series, where the small parameter is the ratio of the surface wave amplitude to the wavelength of the radio wave. Ideally, one should obtain a solution in the form of a complex scattering matrix, but in standard practice, and in view of the stochastic properties of the surface representation, one solves instead for the mean scattered intensity, almost invariably in a linear polarization basis.

In the following section, we review some well-known results from HF scattering theory, pointing out some features that are seldom mentioned explicitly because the usual context—remote sensing of the ice-free ocean surface—does not require us to reflect upon the ‘obvious’ assumptions. As we shall see, these features open the way to our objective, namely, a means of characterizing sea ice from its HF radar signature.

3.2.1. First-Order Theory

At the first-order in the expansion, the expressions for the bistatic scattering coefficients obtained with the small perturbation method take the general form [28]

$${\sigma }_{pq}^{\left(1\right)}({\widehat{k}}_{scat},{\widehat{k}}_{inc})=2^4\pi k_0^4{|G({\widehat{k}}_{scat},{\widehat{k}}_{inc})|}^2{|{\chi }_{pq}|}^{2}S({\overrightarrow{\kappa }}_B)$$

(1)

where $k_0 \equiv |{\overrightarrow{k}}_{inc}| \cong |{\overrightarrow{k}}_{scat}|$ is the wavenumber of the illuminating radio wave, $G\left({\widehat{k}}_{scat},{\widehat{k}}_{inc}\right)$ is a purely geometric function that takes different forms for skywave and surface wave radars, ${\chi }_{pq}$ embodies the polarization-dependence of the scattering amplitude, and $S\left({\overrightarrow{\kappa }}_B\right) \equiv S\left({\overrightarrow{\kappa }}_B\left({\widehat{k}}_{scat},{\widehat{k}}_{inc}\right)\right)$ is the spatial power spectral density of the surface evaluated at the Bragg-resonant wavevector ${\overrightarrow{\kappa }}_B$ given variously by

$$\overrightarrow{\kappa } \equiv {\overrightarrow{\kappa }}_B=\left(1-\widehat{z}\widehat{z}·\right)[{\overrightarrow{k}}_{scat}-{\overrightarrow{k}}_{inc}]$$

(2)

$$\begin{gathered} {\overrightarrow{\kappa }}_B=\left(-k_0\sin {\theta }_{scat}\cos \left({\varphi }_{scat}-{\varphi }_{inc}\right) \\ +k_0\sin {\theta }_{inc} , k_0\sin {\theta }_{scat}\sin \left({\varphi }_{scat}-{\varphi }_{inc}\right) \right) \end{gathered}$$

(3)

The scattering geometry is shown in Figure 3.

Equation (1) can be viewed as the reduction of the integral form

$${\sigma }_{pq}^{\left(1\right)}\left({\widehat{k}}_{scat},{\widehat{k}}_{inc}\right)=2^4\pi k_0^4{|G\left({\widehat{k}}_{scat},{\widehat{k}}_{inc}\right)|}^2{|{\chi }_{pq}|}^2{{\displaystyle \int }}^{}S\left(\overrightarrow{\kappa }\right)\delta \left(\overrightarrow{\kappa }-{\overrightarrow{\kappa }}_B\right)d\overrightarrow{\kappa }$$

(4)

where the delta function implements the Bragg condition.

The scattering coefficients defined by Equation (4) contain no explicit information about the ocean wave dynamics, as manifested in the temporal development of the radar echoes, most commonly expressed in the Doppler domain (although wavelet analysis and other time-frequency decompositions have their niche applications). For space-borne microwave radars, the Doppler domain is seldom exploitable for remote sensing purposes, but for HF radars, Doppler information plays the central role in all geophysical investigations.

The extension of the scattering theory to time-dependent surfaces is accomplished by identifying the components that make up the spatial power spectrum with progressive waves whose spatial and temporal properties are linked by the dispersion relation that embodies the physics of the medium. In the case of small amplitude surface gravity waves on the free ocean surface, the relevant hydrodynamic equation is the linearized inviscid Euler equation, which leads to the familiar dispersion relation

$$\omega =\mathsf{\Omega }\left(\overrightarrow{\kappa }\right)=\sqrt{g|\overrightarrow{\kappa }|\tanh \left(|\overrightarrow{\kappa }|H\right)} \to \sqrt{g|\overrightarrow{\kappa }|} as H\to \infty$$

(5)

Here $g$ denotes the acceleration due to gravity and H is the water depth. For the case of first-order radar scatter, where only the free Bragg-resonant waves contribute, the observed frequency in a geocentric frame is given by

$${\omega }_{obs}=\mathsf{\Omega }\left({\overrightarrow{\kappa }}_B\right)-\overrightarrow{\kappa }·\overrightarrow{U}$$

(6)

where the intrinsic wave frequency $\mathsf{\Omega }\left({\overrightarrow{\kappa }}_B\right)$ is augmented by the term $\overrightarrow{\kappa }·\overrightarrow{U}$, with $\overrightarrow{U}$ denoting the local bulk motion of the water body due to any prevailing current.

Left unstated in this simple but widely adopted formulation is the assumption of barotropic current flow. Once baroclinic effects appear, the shear vorticity renders the potential theory formulation of wave motions invalid, and a more sophisticated treatment of the fluid motions is required. The HF radar community has, almost without exception, adopted a heuristic model that might best be described as mean barotropic equivalent flow, where $\overrightarrow{U}$ is identified with the mean of the baroclinic current taken over a depth range proportional to ${\kappa }^{-1}$.

Accordingly, the radar scattering coefficient for first-order scattering from the time-dependent surface takes the form

$$\begin{gathered} {\sigma }_{pq}^{\left(1\right)}\left({\widehat{k}}_{scat},{\widehat{k}}_{inc},\omega \right)=2^4\pi k_0^4{|G\left({\widehat{k}}_{scat},{\widehat{k}}_{inc}\right)|}^2{|{\chi }_{pq}|}^2\times \\ {\displaystyle \sum }_{m=\pm 1,-1}{{\displaystyle \int }}^{}S\left(\overrightarrow{\kappa },\omega \right) \delta \left(\overrightarrow{\kappa }-m{\overrightarrow{\kappa }}_B\right)· \delta \left(\omega -m\mathsf{\Omega }\left({\overrightarrow{\kappa }}_B\right)+\overrightarrow{\kappa }·\overrightarrow{U}\right)d\overrightarrow{\kappa }d\omega \end{gathered}$$

(7)

where now there is a dispersion relation constraint implemented by another delta function, supplementing the Bragg condition.

Numerous experiments have demonstrated that the selectivity of the Bragg scattering mechanism in HF radar observations of the sea surface is remarkable. Figure 4 shows one such example, recorded with a low power (15 W) HFSWR in Gulf St Vincent, Australia, in 1995.

As in most instances, the spatial resonance is seen here as two peaks corresponding to $m=\pm 1$, that is, one from the waves approaching the radar and one for the waves receding, with positive and negative Doppler shifts, respectively. If the ambient current $\overrightarrow{U}$ presents a non-zero component V aligned with ${\overrightarrow{\kappa }}_B$,

$$V=\overrightarrow{U}·\frac{[{\overrightarrow{k}}_{scat}-{\overrightarrow{k}}_{inc}]}{|{\overrightarrow{k}}_{scat}-{\overrightarrow{k}}_{inc}|}$$

(8)

A displacement of the frequencies of the advancing and receding waves occurs, as indicated in the figure. For this measurement, a CIT of 100 s was employed, and it is apparent that the limit to the spectral width of the resonance peaks is set at least as much by the CIT and the radar waveform bandwidth as by the selectivity of the Bragg scattering mechanism. This selectivity is the key to the standard mission of measuring ocean currents, and accuracy is frequently cited as ±5–8 cm/s. Perhaps it is worth noting that assigning a single numerical value to current speed, or even to measurement accuracy, is, in a sense, misleading, as the surface current velocity, visualized as a point-wise defined quantity, will inevitably vary over a resolution cell in accordance with the turbulence commanded by the Navier–Stokes equation. Nevertheless, it is a most convenient fiction that maps directly onto the intuitive notion of bulk fluid motion.

3.2.2. Second-Order Theory

Valuable as the first-order scattering echoes are, they yield information related to only two waves of the continuum of waves present on the sea surface. The hierarchy of multiple Bragg scattering processes also returns echoes that are observable when the radar has sufficient sensitivity, and these echoes yield information about all the waves on the surface. In the simplest of these processes, Bragg scattering from one wave train directs the incident radio wave not towards the radar receiver but onto a second wave train where another Bragg scattering event occurs. For every instance where the combination of the two Bragg scattering events has the resultant effect of directing the double-bounce scattered field towards the radar receiver, energy will appear in the Doppler spectrum, at a Doppler frequency given by the sum of the intrinsic frequencies of the two participating surface gravity waves,

$$\omega =\mathsf{\Omega }\left({\overrightarrow{\kappa }}_1\right)+\mathsf{\Omega }\left({\overrightarrow{\kappa }}_2\right)$$

(9)

while the sum of the wavevectors of the two participating surface gravity waves must satisfy

$$\overrightarrow{\mathsf{{\rm K}}}\equiv {\overrightarrow{\kappa }}_1+{\overrightarrow{\kappa }}_2=\left(1-\widehat{z}\widehat{z}·\right)[{\overrightarrow{k}}_{scat}-{\overrightarrow{k}}_{inc}]$$

(10)

where $\widehat{z}$ is the unit normal to the surface.

Having extended the electromagnetics model to the second order, we must do the same for the hydrodynamics. At the second order, a propagating disturbance associated with a wavevector $\overrightarrow{\mathsf{{\rm K}}}$can be expressed as an integral over pairs of first-order waves. In this case, the Bragg condition requires that

$$\overrightarrow{\mathsf{{\rm K}}}=m_1{\overrightarrow{\kappa }}_1+m_2{\overrightarrow{\kappa }}_2$$

(11)

with $m_1, m_2=\pm 1$as above, while the dispersion relation constraint requires that

$$\mathsf{\Omega }(\overrightarrow{\mathsf{{\rm K}}})=m_1\mathsf{\Omega }({\overrightarrow{\kappa }}_1)+m_2\mathsf{\Omega }({\overrightarrow{\kappa }}_2)$$

(12)

with a physics-dependent coupling coefficient. However, a problem arises with waves on the free ocean surface as the dispersion relation (5) does not support free waves satisfying both (11) and (12). Instead, the resultant disturbance must exist as an evanescent wave, unable to propagate independently but phase-locked onto its parent waves, where it can still generate its own strong radar signature. This raises an interesting question: although triad interactions with energy exchange are not permitted on free ocean surfaces, are they allowed on ice-covered waters? We shall return to this question later.

Combining the second-order effects, a Doppler continuum is generated, with its shape dependent on the amplitudes of all the gravity waves on the surface. The solution for the resultant second-order Doppler spectrum was developed in [29,30,31,32,33]; it takes the general form

$$\begin{gathered} {\sigma }_{pq}^{\left(2\right)}({\overrightarrow{k}}_{scat},{\overrightarrow{k}}_{inc},\omega )=2^4\pi k_0^4{|G|}^2{|{\chi }_{pq}|}^2\times \\ {\displaystyle \sum }_{m_1,m_{2=+1,-1}}{{\displaystyle \iint }}^{}d{\overrightarrow{\kappa }}_{1}d{\overrightarrow{\kappa }}_{2}S(m_1{\overrightarrow{\kappa }}_1)S(m_2{\overrightarrow{\kappa }}_2)\Gamma ({\overrightarrow{k}}_{scat},{\overrightarrow{k}}_{inc},m_1{\overrightarrow{\kappa }}_1,m_2{\overrightarrow{\kappa }}_2,\omega )\times \\ \delta ({\overrightarrow{\kappa }}_1+{\overrightarrow{\kappa }}_2-\left(1-\widehat{z}\widehat{z}·\right)[{\overrightarrow{k}}_{scat}-{\overrightarrow{k}}_{inc}]). \delta (\omega -m_1\mathsf{\Omega }({\overrightarrow{\kappa }}_1)-m_2\mathsf{\Omega }({\overrightarrow{\kappa }}_2)-\overrightarrow{K}·\overrightarrow{U}) \end{gathered}$$

(13)

where the kernel $\Gamma ({\overrightarrow{k}}_{scat},{\overrightarrow{k}}_{inc},m_1{\overrightarrow{\kappa }}_1,m_2{\overrightarrow{\kappa }}_2,\omega )$ embraces both electromagnetic and hydrodynamic contributions. Expressions for $\Gamma$ as applicable to ice-free ocean surfaces take the form

$$\Gamma ={|{\Gamma }_{EM}+{\Gamma }_{HYD}|}^2$$

(14)

where

$${\Gamma }_{EM}=\frac{[\left(1-\widehat{z}\widehat{z}·\right){\overrightarrow{k}}_{inc}]·[\left(1-\widehat{z}\widehat{z}·\right){\overrightarrow{k}}_{scat}]·[\left(1-\widehat{z}\widehat{z}·\right){\overrightarrow{k}}_{inc}-k_0]}{2k_0{\cos }^2\left({\phi }_{scat}-{\phi }_{inc}\right)\left(\sqrt{[\left(1-\widehat{z}\widehat{z}·\right){\overrightarrow{k}}_{inc}]}·\left[[\left(1-\widehat{z}\widehat{z}·\right){\overrightarrow{k}}_{inc}]-2k_0\right]+i{k}_0\Delta \right)}$$

(15)

$$\begin{gathered} {\Gamma }_{HYD}=-\frac{i}{2}\left[\frac{\left({\overrightarrow{\kappa }}_1·{\overrightarrow{\kappa }}_2-{\kappa }_1{\kappa }_2\right)\left({\omega }^2+{\omega }_B^2\right)}{m_1{m}_2\sqrt{{\kappa }_1{\kappa }_2 }\left({\omega }^2-{\omega }_B^2\right)}+{\kappa }_1+{\kappa }_2 \\ +\frac{\omega \left[{\left(m_1\sqrt{g{\kappa }_1}\right)}^{3}csc{h}^2\left({\kappa }_{1}H\right)+{\left(m_2\sqrt{g{\kappa }_2}\right)}^{3}csc{h}^2\left({\kappa }_{2}H\right)\right]}{g\left({\omega }^2-{\omega }_B^2\right)}\right] \end{gathered}$$

(16)

with H as the water depth [34] and $\Delta$ the normalized surface impedance given by

$$\Delta =\frac{1}{N}\sqrt{1-\frac{1}{N^2}}\mathrm{with}N=\sqrt{{ϵ}_r+i\frac{\sigma }{{ϵ}_0{\omega }_{radar}}}$$

Here ${ϵ}_r$ is the relative permittivity of sea water, σ its conductivity, ${ϵ}_0$ the permittivity of free space and ${\omega }_{radar}$ the angular frequency of the radio wave.

3.3. Remote Sensing of the Ocean Surface with HF Radar

3.3.1. Free Ocean Surface

The solution of the forward problem of radio wave scattering from the sea surface, as expressed in Equations (7) and (13), is of interest primarily because it leads to the inverse problem of extracting the sea surface parameters from radar echoes. The most accessible parameter, retrievable from even low power radar echoes, is the ${\overrightarrow{\kappa }}_B$-aligned component of the surface current vector. Deploying two radars to illuminate the same patch of sea from different directions yields the full velocity vector $\overrightarrow{U}$, and this has turned out to be by far the most successful, widespread, and robust application of HF surface wave radars. Inversion of the Doppler spectrum including both first- and second-order echoes to extract $S(\overrightarrow{\kappa })$ is a more challenging task, especially as it requires data of a high dynamic range, which are not always available with low power radars. Techniques to achieve this fall into several categories, as outlined in [35]. Further, to the extent that $S(\overrightarrow{\kappa })$ constitutes a full local description of the ocean surface, we can relax the data requirements somewhat to retrieve integral quantities via the extraction of robust features of the Doppler spectrum, in conjunction with empirical relationships. Specifically, writing $S(\overrightarrow{\kappa })$ in polar form, we identify the hierarchy

$$\begin{array}{ccc}S\left(\kappa ,\phi \right)& \Rightarrow S\left(\kappa \right)& ={{\displaystyle \int }}^{}S\left(\kappa ,\phi \right)J\left(\kappa ,\phi \right)d\phi \\ & \Rightarrow H_s& =4{{\displaystyle \iint }}^{}S\left(\kappa ,\phi \right)J\left(\kappa ,\phi \right)d\phi d\kappa \end{array}$$

(17)

projecting the directional spectrum onto the nondirectional spectrum $S\left(\kappa \right)$ and thence to the significant waveheight, $H_s\equiv 4\sqrt{{\eta }^2}$ where $\eta \equiv \eta (\overrightarrow{r})$ is the surface displacement and the kernel $J\left(\kappa ,\phi \right)$ embodies the appropriate Jacobian.

These procedures have been implemented in many HF radar systems. Techniques for modelling and retrieving some more exotic surface phenomena, such as the presence of surfactants, wave damping due to rainfall, and ship wakes, have been developed but these are not in general use. On the other hand, detection of echoes from ships, once associated only with military HF radars, has become a staple of commercial remote sensing radars.

3.3.2. Extension to Ice-Covered Surfaces

When we turn to the situation where ice is present, it seems logical to adapt and apply the techniques that have proved so successful for ice-free surfaces. Accordingly, we start with the premise that an analysis of measured Doppler spectra resulting from motions of the ice-covered surface, as it responds to disturbances propagating through the ice field, should enable us to identify the type of ice and, ideally, measure its parameters. As we shall see, this expectation is met, at least in part, though the limits of the approach, and the demands it places on the data have not yet been fully mapped. In the following section, we shall examine the types of ice that may be encountered and identify the dispersion relations of wave-like disturbances in these media. Following [36], we shall refer to these disturbances as ice-coupled waves. This terminology is adopted to remind ourselves that, while the waves of interest originate as free surface gravity waves in the open sea, they are qualitatively different when they enter the ice-water environment.

Before we review the physics of ice-coupled waves, derive the dispersion relations, and substitute them into the scattering formulae, we need to check that other assumptions built into the free surface remote sensing methodology are not violated when ice is present. One question in particular needs to be addressed: is the small perturbation method for HF radio wave scattering still valid?

The primary criterion of validity is a constraint on the roughness scale, $k_0\eta \ll 1$. In the case of the free surface, there is no ambiguity: the water surface displacement under the wave action defines the roughness scale. For ice-covered surfaces, we are still concerned with the displacement due to ice-coupled waves, but there are additional forms of surface roughness present, associated with (i) the ice surface, and (ii) the underside of the ice, which coincides with the upper surface of the sea water medium. These introduce some electromagnetic subtleties associated with the layered structure, including the possibility of additional electromagnetic wave modes, but the key point to remember is that the Doppler spectrum results from the temporal modulation of the phase of the radio wave, and this arises solely from the time-varying geometry associated with the ice-coupled wave. The intrinsic ice roughness, observable as the ice freeboard and the ice thickness, along with its electrical permittivity and conductivity, will manifest itself almost exclusively in the spatial characteristics of the scattered field, and hence the echo intensity, as long as the ice upper surface motion closely follows that of the local ice-coupled wave as defined by the ice lower surface. The relative magnitudes of the intrinsic and ice-coupled wave roughness scales will depend entirely on the circumstances, but for illustration it is useful to consider Figure 5, which shows an Envisat altimeter profile of the sea ice freeboard over a path of 2500 km, adapted from [37]. Seldom does the freeboard reach 0.3 m; the mean over the entire record is 0.162 m. This is consistent with typical visual images of first-year Arctic sea ice. Note, though, that with ice thickness roughly nine times freeboard, there is a substantial variation in ice thickness and hence ice hazard for human activities, such as ship transit.

From an informal survey of other reports, the freeboard of one-year ice generally lies in the range 0–0.2 m [38], while for multi-year ice, it falls in the range 0.1–0.3 m [37]. (Ice thickness is, of course, an order of magnitude greater.) For HF radio waves, the wavelength varies from 100 m at 3 MHz to 10 m at 30 MHz, though we note that preferred frequencies lie between 3 MHz and 10 MHz. Thus, $k_0\eta$ as defined by the upper surface of the ice typically takes values in the range 0.0–0.05. As for ice-coupled wave amplitudes, these can be substantial; [39] reported waves of 3 m amplitude at a location more than 350 km inside the ice field and other reports stretch to 500 km.

A related consideration, often overlooked, is an implicit assumption regarding the surface slope, arising from the use of the Taylor expansion to transfer the boundary conditions to the mean surface. Clearly there are sharp gradients at the edges of ice floes, though not in the physical models adopted for thin elastic sheets and viscous layers, but the convergence of the perturbation expansion involves the product of the gradient and the displacement, so for the very small displacements associated with almost all parts of the ice-covered surface, some tolerance of local prominences in slope can be accepted.

We must also examine the extent to which the kernels of the scattering integrals change when ice is present. It seems reasonable to assume that the electromagnetic term ${\Gamma }_{EM}$ developed in the original Rice formulation will apply, with a suitable choice for the electrical permittivity, but we cannot expect the nonlinear interaction of ice-coupled waves to take exactly the same form as that for free surface waves, except for the limiting case of a vanishingly small ice cover. A model for this nonlinear interaction term is presently under development so, for the present, we shall use the free surface hydrodynamic kernel. Although the resulting Doppler spectra may differ somewhat from those computed with the correct kernel, the object of our study is the dependence on ice species and ice parameter values, and these dependences should reveal themselves irrespective of slight variations in the hydrodynamic kernel. Indeed, it must be remembered that the first-order echoes, which are so informative, do not rely on these kernels.

Finally, another (somewhat speculative) prospect that contrasts with the free surface situation is that the medium may be anisotropic for ice-coupled waves. During the break-up phase, the fracturing of ice sheets and large floes into smaller floes seems to occur preferentially during the passage of crests and (perhaps) troughs of propagating ice-coupled waves, which are usually highly directional. The physical explanation for this phenomenon lies in the observation that the fractures occur along the loci of maximum surface curvature, where flexural stress peaks, so a directional wave spectrum incident on the ice field would engender this anisotropy.

Normally, the floe field is represented by a scalar floe size distribution. Visual imagery of such ice fields frequently shows remarkably polygonal floe shapes, suggesting that additional parameters are needed to describe the distribution. One attribute that, to our knowledge, has escaped detailed analysis is the edge orientation distribution. We would expect this to impact not only on the properties of floe–floe collisions and wave–floe scattering, but it would also change the electromagnetic properties of the surface, as the network of open leads makes the surface electrically analogous to a wire grid on a substrate.

In order to assess the anisotropy of ice field geometry, we have developed tools to automate the extraction of floe geometry from aerial photographs. The methodology is based on the empirical observation that the shape of most floes can be well approximated by polygons of an order of less than or equal to eight, affording a computationally efficient means of evaluating the distribution of edge orientations. Figure 6 presents a typical input-output result of a field of ‘polyfloes’; we have yet to complete the development of tools for analyzing and characterizing anisotropy so we shall not pursue this matter here.

On the basis of these assessments, we conclude that the perturbation theoretic expansion for the scattered field can be applied to compute the Doppler spectrum of ice-coupled waves once the appropriate dispersion relations have been imposed.

3.4. Previous Studies and Deployments of HF Radars in the Ice Zone

To answer the second question raised at the beginning of this section, investigations into the possible utility of HF radars for ice monitoring date back to the 1970s, when skywave radars looking northward from the continental US observed reflections and forward scattering from the polar ice cap. These studies were performed mainly under the aegis of the Polar Fox research program [40], using a facility in Maine with the flexibility to conduct measurements across almost the entire HF band. The primary objective was not to study the geophysical properties of the Earth’s surface but to assess the utility of these ground echoes for calibrating the radar sensitivity. Experimental measurements yielded estimates of the static scattering coefficients of (i) open sea, (ii) sea ice, and (iii) the Greenland ice cap, which researchers at the time interpreted as being adequate for calibration [41].

During the early 1980s, Canadian interest in the potential applications of HF radar at high latitudes developed, primarily at Memorial University, leading to the deployment of a US CODAR radar in Labrador, in 1984. This system successfully detected icebergs and validated a theoretical model of iceberg radar cross section developed by researchers at Memorial University [42,43].With commercial and government support, this initiative led to the deployment of a large HFSWR system at Cape Race on the coast of Newfoundland; it became operational in 1990 and extended the range at which icebergs were detected to over 200 km. The transmit and receive arrays are pictured in Figure 7 and Figure 8. The radar operated at frequencies in the band 3–5 MHz.

In parallel with this program, but also involving Memorial University, industry, and government, a second radar was constructed in Cape Bonavista, Newfoundland, in 1990. Operating at 1.95 MHz, which lies in the MF band, just below HF, icebergs were detected and tracked at ranges nearing 300 km and extrapolation suggested the possibility of the detection of very large icebergs at a range of 500 km. Theoretical modelling and assessment of the radar data indicated that each of these radars was able, in principle, to detect the edge of a solid ice field. Based on our understanding of these radars and the theory presented in this paper, we can assert that, in principle, each could be used for the detection of ice-coupled waves if located near the marginal ice zone. Figure 9 and Figure 10 show the Cape Bonavista transmitting and receiving antenna arrays.

The expansion of offshore resource extraction in the Beaufort Sea in the 1980s led to the deployment of a CODAR radar on Cross Island in 1984. This system measured ocean currents and tracked transponders mounted on small floes, and in addition observed the breakup of a large ice pack, initially almost stationary but accelerating as it broke up under the action of thermal and wind stresses, opening up regions of exposed sea through which the small floes moved with the same velocity as the local current. These measurements were made at a range of 15 km, consistent with the very low power of the CODAR and the attenuation of the relatively high HF frequency (25.4 MHz) over ice. Since that time, CODAR systems have been deployed by US researchers in numerous other Arctic and Antarctic locations, focusing on the measurement of ocean currents [44,45,46,47,48,49]. The low power and low spatial resolution of these radars reduces their utility for the ice-coupled wave monitoring mission, but their ultra-compact design and high Doppler resolutions—able to estimate velocities to ±2–4 cm/s—count in their favor. Figure 11 shows a CODAR SeaSonde near Prudhoe Bay, Alaska, with Figure 12 showing a SeaSonde radar at Anvers Island, Antarctica. Perhaps the earliest current measurements with HF radar in the Antarctic are those reported in [50], using a UK OSCR-II radar. None of these investigations attempted to study the echoes from sea ice with a view to the kind of analysis we introduce in the present paper.

The question of the reliability of ice edge detection with HF radar was addressed by Japanese researchers in 2013, using data from 2004 which overlapped with C-band radar measurements [51]. Using a CODAR system, they mapped the ice edge north of Hokkaido where it varies rapidly under strong currents associated with the Soya warm current. In their approach, it was simply the absence of radar echo that determined the onset of the ice-covered sea, exploiting the higher attenuation of high HF band signals over ice. Their results affirmed the performance of HF radar in this role, to a range of over 30 km when the ice field is viewed from the open sea side.

A very different experiment was conducted by Canadian researchers in the lower St Lawrence Estuary, Quebec [52,53]. They set out to compare the performance of two HFSWR systems—one CODAR and one German WERA—in the mission of measuring currents in partially ice-covered waters. We might remark, incidentally, that the results show similar dependence on ice concentration, but the noteworthy feature of this work from our perspective is the attention that was paid to the stages of development of the ice. Pictures of the phased array antennas of the WERA installation are presented in Figure 13 and Figure 14. The latter confirms the feasibility of mounting antennas directly onto the ice sheet.

Finally, a recent paper by an international team extended the earlier study of the Soya current, re-examining the 2004 datasets with different processing techniques and validating the HF radar measurements of ice drift velocity by comparisons with ADCP and drifting buoy data [54]. Again, the physical properties of the ice played no role.

It is evident that many researchers have recognized the potential of HF radar to provide operationally useful information on the spatial distribution of sea ice and its drift under the action of ocean currents. Decades of high-quality research speak to this. Yet, and this comes as something of a surprise, until the investigation reported here, ice-related observations and theoretical studies have dealt only with detecting the ice edge, detecting icebergs, and measuring ocean currents where ice may be present. The idea of measuring the dynamical response of sea ice appears to have been overlooked, or more likely, deemed to lie beyond prevailing radar capabilities. In the remainder of this paper, we shall develop the theoretical framework within which this objective might be achieved.

4. Sea Ice: Forms and Properties

4.1. Formation and Evolution of Sea Ice

When the air temperature drops below about −1.8 °C, depending on salinity, minute ice crystals form in the surface layer of the sea. As these frazil particles grow, the water becomes slushy, and then, as the particles coalesce, the viscosity increases sharply, forming grease ice. What happens next depends on the level of wave activity. In calm conditions the surface eventually freezes entirely, forming nilas, a thin layer of sheet ice. This provides the base for coagulation, where the ice thickens, initially rather quickly but then more slowly as the ice layer itself impedes the heat flux from water to air. This layer is termed young ice when the thickness lies in the range 10–30 cm. Over one winter, young ice may grow to as much as 2 m in thickness. This first-year ice is of highest importance to human activities. At any stage, the young ice field may be subjected to wave action from waves generated outside the ice zone and propagating through it. This results in the break-up of the ice field into floes, formally defined as having a dimension greater than 20 m, which may then drift under prevailing currents or wind forcing.

Under the other possibility, where the grease ice evolves in the presence of wave action, the frozen material organizes itself into centimeter-scale clumps called shuga, that collide and fuse into pancake ice, which are flat, disc-shaped structures. Under further wave action, the pancakes may be subject to collisions, rafting and subsequent conglomeration into floes. These pathways are illustrated in Figure 15.

Mature ice fields can develop a complex topography in the form of pressure ridges, leads, polynyas, melt ponds, bergs, and other features, but for now we shall assume homogeneous surfaces satisfying the SPM, so that the notion of a surface perturbed by propagating ice-coupled waves is appropriate.

4.2. Properties of Sea Ice

From the perspective of HF radar, the key attributes of interest are (i) the electrical properties that govern the electromagnetic scattering process, and (ii) the mechanical and structural properties that determine the spatial and temporal response of the surface to waves progressing through the ice-ocean system. Together, these guide the development of models of the ice that can be used to explore its dynamical behavior on short time scales and to compute the radar signature. Processes such as ice formation and melting take place on longer time scales, where they are controlled by thermodynamic properties.

4.2.1. Electrical Properties

The electrical properties have essentially no bearing on the fluid characteristics of ice-coupled waves, but they have a strong bearing on the interaction of incident radio waves with the ice field. Accordingly, the permittivity, permeability and conductivity of sea ice have been measured by many researchers over the past century, though with more emphasis on frequencies above and below the HF band. During our investigation, the data in [55,56,57] have proven useful. In addition to mapping generic behavior, such as the temperature dependence of conductivity, one of the most significant observations has been the recognition of the long-term change in electrical properties as the small pockets of brine trapped in the newly-formed ice gradually leach away, thereby removing the most highly conducting component. This fact has been exploited by SAR and other microwave sensors to distinguish first-year ice from multi-year ice. A second complication arises from the fact that the brine pockets are not uniformly distributed in their orientation, resulting in the need to employ a dielectric tensor to describe the interaction with electromagnetic waves. In the present context, it is the effect of the permittivity and conductivity on HF surface wave propagation that has the most dramatic operational impact, though the electrical parameters also act through the polarization terms in the scattering model equations. In Figure 16, we have plotted the ITU standard values of the permittivity and conductivity of sea ice, together with indications of the disparities with the corresponding values for seawater.

4.2.2. Mechanical and Structural Properties

It scarcely needs to be stated that the diverse physical structures assumed by sea ice demand their own phenomenological models, variously emphasizing macroscopic or microscopic features, solid vs. liquid properties, lossy or loss-free distortion, linear or nonlinear behavior, local or nonlocal aspects, and the identification of dominant modes when a multiplicity of wave-like responses can co-exist [58,59,60].

From the dynamical perspective that is of concern to this study, we are interested in the spectrum of the wave-like excitations of the composite system, which consists of an ice layer floating on sea water. The ice may take any of the forms discussed above, or a mixture of these, so, in principle, a mathematical model describing the excitations must be able to account for a diverse set of mechanical and rheological properties. In particular, researchers have identified the roles of inertial mass, elasticity (in tension, compression and shear), flexural rigidity, viscosity, porosity and fracture toughness.

While there has been considerable progress recently in the formulation of models that are able to address multiple attributes within a physically and mathematically consistent framework, it has been observed in practice that the dynamics of the most common forms of sea ice are quite well represented by models involving just one or two dominant constitutive parameters. This has the practical advantage of facilitating inversion of the observational data to retrieve the ice parameters. A priori information in the form of focused measurements that place bounds on parameter ranges also serve to assist retrieval [61,62,63].

The best-known ice model classes, which we might call the canonical models, are:

• The mass loading model;
• The thin elastic plate model;
• The viscous layer models;
• The viscoelastic layer models;
• The discrete floe models.

All of these except the last can be modelled by linear partial differential equations with constant coefficients on an unbounded domain. Such equations admit plane wave solutions of the form $\eta \left(x,t\right)=e^{i\left(\kappa x-\omega t\right)}$ with the constraint that the values of $\kappa$ and $\omega$ must satisfy a dispersion relation, $\omega =\mathsf{\Omega }\left(\kappa \right)$ such as the one discussed in Section 3.2. This means that almost all of the ice models listed are amenable to representation in our scattering theory to yield a Doppler spectral signature.

The remaining case—a distribution of discrete ice floes of dimensions comparable with, or greater than, the wavelength of the penetrating wave—poses many problems. Does it make sense to define a dispersion relation in this case? As an HF radar resolution cell will contain many floes, one is tempted to believe that the Fourier transform of the surface displacement field $\eta \left(x,t\right)$ into the $(\omega ,\overrightarrow{\kappa })$ spectral domain will reveal a concentration of energy along some one-dimensional curve, just as it does for the other model classes, but this is by no means guaranteed. We must also take into account that the individual floes have their own intrinsic dynamical responses to excitation, which one would expect to be functions of the floe parameters, so the notion of a homogeneous ice field with invariant response characteristics must be abandoned at some length scale. If the floes occupy only a small fraction of the surface, perhaps their impact on the $\left(\omega ,\kappa \right)$-space distribution of energy might be merely to modify the angular spectrum of waves propagating on the ambient fluid surface waves, presumably described by a viscous layer model or even a free surface model.

Our goal in this study is limited to showing how different models predict different HF radar signatures and thence to establish whether the differences might be sufficient to yield an operational ice diagnosis and measurement capability. As long as our chosen models faithfully embody the main dynamical characteristics of the corresponding ice forms, we will be able to achieve our objective, but, clearly, we are unable to do this for the discrete floe models. Even with the other models, for which we do possess the dispersion relations, we must keep in mind that, apart from their intrinsic inadequacies, these models represent ‘pure states’ that do not span the space of possible sea ice conditions as mixtures of ice types may exist within a radar resolution cell.

4.3. Dispersion Relations for Ice-Coupled Waves

Two distinct approaches to modelling ice-coupled waves have been explored by the ice modelling community. The first adopts a single fluid–atmosphere boundary model, just as in the free surface case, and solves for progressive waves, modifying the dynamic boundary condition at the surface to simulate the stresses imposed by the prescribed ice rheology. While simplistic in form, models in this category have yielded satisfactory results in some ice scenarios and have the advantages of explicit form and, in some cases, invertibility. In the second approach, a finite thickness ice layer is explicitly represented, so two interfaces are involved; the ice and the underlying fluid are endowed with their respective constitutive properties and solutions are sought in the form of progressive waves. We shall not include detailed derivations of all the principal models here—for those we provide ample references—but it greatly adds to one’s understanding of the limitations of the models if one steps through an outline of the theory behind the simplest classes, seeing where the approximations are made. It may even aid the development of customized generalizations, such as the incorporation of a snow layer. Rather than relegate this outline to an appendix, we address it here.

Accordingly, in the first approach, the mathematical description of ice-coupled waves begins with the specification of the following:

(i) The Navier–Stokes equation to model the fluid motion;

(ii) Kinematic boundary conditions for the interfaces;

(iii) Dynamic boundary conditions for the interfaces;

while in the second, these must be supplemented with

(iv) Equations to model the stress-strain tensors of the ice medium.

The Navier–Stokes equation can be written in the general form

$$\rho \left(\frac{\partial \overrightarrow{v}}{\partial t}+\overrightarrow{v}·\nabla \overrightarrow{v}\right)=-\nabla p+\nabla ·\left[\mu \left(\nabla \overrightarrow{v}+{\left(\nabla \overrightarrow{v}\right)}^T-\frac{2}{3}\left(\nabla ·\overrightarrow{v}\right)\stackrel{=}{I}\right)+\zeta \left(\nabla ·\overrightarrow{v}\right)\stackrel{=}{I} \right]+\rho \overrightarrow{g}$$

(18)

where $\overrightarrow{v}$ is the fluid velocity, $p$ is the pressure, $\mu$ and $\zeta$ are the first and second viscosities, $\stackrel{=}{I}$ is the identity tensor, $\rho$ is the fluid density, $\overrightarrow{g}$ is the gravitational acceleration and $\overrightarrow{M}$ represents any other body forces that might be present. (It might be thought that it should suffice to employ the simplified form $\mu {\nabla }^2\overrightarrow{v}$ instead of the complicated viscosity term in the square brackets, or else to assume inviscid flow, but the frazil ice crystals that develop into grease ice fields tend to form sometimes as needles and often as thin discs; in each case the viscosity of an assembly of such structures presents an anisotropic resistance and hence a proper treatment requires a tensor representation, as reported in [64].) All that said, we can establish the essential characteristics of the simpler canonical models of ice-coupled waves in an elementary but informative way. First, assume the flow is inviscid and irrotational, so we can represent the velocity by a scalar potential,

$$\overrightarrow{v}=\nabla \phi$$

(19)

Then, further assuming that the fluid is incompressible, $\phi$ must satisfy the Laplace equation,

$${\nabla }^2\phi =0$$

(20)

In the absence of viscosity, the Navier–Stokes equation reduces to a Euler-type equation,

$$\rho \left(\frac{\partial \overrightarrow{v}}{\partial t}+\overrightarrow{v}·\nabla \overrightarrow{v}\right)=-\nabla p+\rho \overrightarrow{g}$$

(21)

Now, integrate once with respect to the spatial coordinates, yielding the Bernouilli equation,

$$\frac{\partial \phi }{\partial t}+\frac{1}{2}{\left(\nabla \phi \right)}^2=\frac{p}{\rho }+g\eta$$

(22)

where $\eta$ is the surface displacement. We linearize by retaining only the first-order terms and evaluating at the still water surface $\eta =0$, yielding an expression for the pressure in the fluid at the interface with the ice,

$$p_{int}=-\rho \left(\frac{\partial \phi }{\partial t}+g\eta \right)$$

(23)

At the same location, the ice dynamics are governed by its equation of motion,

$$p_{int }-\tilde{M}\eta ={\rho }_{ice}h \frac{{\partial }^2\eta }{\partial t^2}$$

(24)

where $\tilde{M}$ is an operator that embodies all the stresses and body forces within the ice that oppose the deformation, ${\rho }_{ice}$ is the density of ice and $h$ is the ice thickness. Equating the two expressions for the pressure, we obtain the dynamical boundary condition that must be satisfied at the interface,

$$\tilde{M}\eta +{\rho }_{ice}h \frac{{\partial }^2\eta }{\partial t^2}=-\rho \left(\frac{\partial \phi }{\partial t}+g\eta \right)$$

(25)

We can rearrange this to separate the terms that yield powers of $\omega$ when applied to harmonic waves, from those that yield powers of $\kappa$,

$${\rho }_{ice}h \frac{{\partial }^2\eta }{\partial t^2}+\rho \frac{\partial \phi }{\partial t} =-\left(\tilde{M}\eta +\rho g\eta \right)$$

(26)

Within a linearized framework, the ice layer stresses, and body forces take the general form

$${\tilde{M}}_j=Q_j{h}^{m_j}\frac{{\partial }^{n_j}}{\partial x^{n_j}}$$

(27)

The dominant stress for the thin plate model is flexural stress, $L$ for which $n_j=4$ and $m_j=3$ while the compressional stress $K$ has the values $n_j=2$ and $m_j=1$.

To obtain the dispersion relation, consider the case of a single sinusoidal ice-coupled wave,

$$\eta \equiv \eta \left(x,t\right)=a \cos \left(\kappa x-\omega t\right)\equiv a\Re \left(e^{i\left(\kappa x-\omega t\right)}\right)$$

(28)

We first need to find $\phi \left(x,t\right)$ that satisfies the new dynamic boundary condition with this form for $\eta .$From the separation of variables, we obtain the general solution to Laplace‘s equation,

$$\phi =\left[A{e}^{\kappa z}+B{e}^{-\kappa z}\right]e^{i\left(\kappa x-\omega t\right)}$$

(29)

where the coefficients need to satisfy the boundary conditions. The no-penetration condition at the seabed requires that

$${\left[\frac{\partial \phi }{\partial z}\right]}_{z=-H}=0=\kappa \left(A{e}^{\kappa H}+B{e}^{-\kappa H}\right)e^{i\left(\kappa x-\omega t\right)}$$

(30)

whence $B=A{e}^{-2\kappa H}$, so

$$\phi =A{e}^{-\kappa H}\left[e^{\kappa \left(z+H\right)}+e^{-\kappa \left(z+H\right)}\right]e^{i\left(\kappa x-\omega t\right)} \equiv C \cosh \left(\kappa \left[z+H\right]\right)e^{i\left(\kappa x-\omega t\right)}$$

(31)

Substituting in the dynamic boundary condition and equating to the definition of $\eta$,

$$ag{e}^{i\left(\kappa x-\omega t\right)}=-\left[\frac{{\sum }_j{Q}_j{\left(i\kappa \right)}^j{h}^{m_j}-{\rho }_{ice}h {\omega }^2}{\rho }a{e}^{i\left(\kappa x-\omega t\right)}-i\omega C \cosh \left(\kappa \left[z+H\right]\right)e^{i\left(\kappa x-\omega t\right)}\right]$$

(32)

Rearranging, we solve for $C$ at $z=0$

$$C=\frac{a}{-i\omega \cosh \left(\kappa H\right)}\left[g+\frac{{\sum }_j{Q}_j{\left(i\kappa \right)}^j{h}^{m_j}-{\rho }_{ice}h {\omega }^2}{\rho }\right]$$

(33)

whence

$$\phi =a \frac{{\sum }_j{Q}_j{\left(i\kappa \right)}^j{h}^{m_j}-{\rho }_{ice}h {\omega }^2+g\rho }{-i\omega \rho \cosh \left(\kappa H\right)}\cosh \left(\kappa \left[z+H\right]\right)e^{i\left(\kappa x-\omega t\right)}$$

(34)

The normal fluid velocity at the boundary is defined by the kinematic boundary condition,

$$\frac{\partial \phi }{\partial z}=\frac{d\eta }{dt}=\frac{\partial \eta }{\partial t}+{\nabla }_h\phi ·{\nabla }_h\eta \approx \frac{\partial \eta }{\partial t}$$

(35)

where again we have linearized the equation. Substituting for $\phi$ and $\eta$

$$-a\kappa \sinh \left(\kappa \left[z+H\right]\right)\frac{{\sum }_j{Q}_j{\left(i\kappa \right)}^j{h}^{m_j}-{\rho }_{ice}h {\omega }^2+g\rho }{-i\omega \rho \cosh \left(\kappa H\right)}{e}^{i\left(\kappa x-\omega t\right)}=-i\omega a{e}^{i\left(\kappa x-\omega t\right)}$$

(36)

and evaluating at the mean interface $z=0$ yields the general form of the dispersion relation

$$\frac{1}{\rho }{\displaystyle \sum }_j{Q}_j{\left(i\kappa \right)}^j{h}^{m_j}-\frac{{\rho }_{ice}}{\rho }h {\omega }^2+g=\frac{{\omega }^2}{\kappa }\coth \left(\kappa H\right)$$

(37)

Thus, for a given set of physical stress mechanisms that we assign to the ice, we will obtain a corresponding dispersion relation, as the following examples demonstrate.

In almost all cases, the form of the dispersion relations allows for multiple wave species to exist, which might seem to render our goal unachievable, and even the application of simple models pointless, but on most occasions only one or two solutions lie close to the real axis in complex $\overrightarrow{\kappa }$ space, so the other roots are attenuated rapidly and have little bearing on the time-varying geometrical dynamics that yields the radar signature. (Whether the same is true for nonlinear wave triads is yet to be determined.).

4.3.1. Mass Loading Model

In the mass loading model [65,66,67], we regard the ice layer as a collection of non-interacting point masses, so $Q_j=0 \forall j$ and the dispersion relation is

$${\omega }^2=\frac{g\kappa }{\left( \coth \left[\kappa H\right]+\frac{{\rho }_{ice}}{\rho }h\kappa \right)}$$

(38)

where the effective mass load (per unit area) is expressed by ${\rho }_{ice}h$. In order for the mechanical stresses to vanish, the ice must be in pancakes or small floes whose dimension is much smaller than the wavelength of the ice-coupled wave. If we let the ice thickness go to zero, we regain the free surface dispersion relation. Figure 17 plots the dispersion relation for the mass loading model with equivalent thicknesses of 0.5 m and 1.0 m, together with that for the free ocean surface. Although we are not yet convinced of the need to accommodate anisotropy, here we plot the dispersion surfaces for the two-dimensional case.

4.3.2. Thin Elastic Plate Model

In the simplest case of an ice sheet with integrity, where stresses lead to strains, the property of flexural elasticity is assigned to the ice layer [68,69]. The associated bending force term is proportional to the fourth derivative of the displacement, which leads to the dispersion relation

$${\omega }^2=\frac{g\kappa +\frac{1}{\rho }L{h}^3{\kappa }^5}{\left(\coth \left[\kappa H\right]+\frac{{\rho }_{ice}}{\rho }\kappa h\right)}$$

(39)

For an isotropic material, the modulus of flexural rigidity $L$ is related to the Young’s modulus $E$ and the Poisson ratio ${\nu }_P$ through

$$L=\frac{E}{12\left(1-{\nu }_P{}^2\right)}$$

(40)

We see an interesting property of this dispersion relation: increasing $h$ changes both numerator and denominator in the same direction, so, depending on the relative importance of flexural rigidity and mass loading, the dispersion relation may result in a phase velocity greater or less than that for a free surface. Waves incident on the ice edge would then refract towards the normal or away from it, resulting in different directional spectra and thus different radar signatures.

If the ice thickness $h$ increases, we might expect elasticity in compression to play a role, in which case the dispersion relation is modified yet again [70], taking the form

$${\omega }^2=\frac{g\kappa +\frac{1}{\rho }\left[L{h}^3{\kappa }^5-Kh{\kappa }^3\right]}{\left(\coth \left[\kappa H\right]+\frac{{\rho }_{ice}}{\rho }\kappa h\right)}$$

(41)

with $K$ the compressional elasticity. In practice, the compressibility becomes significant only when the ice layer thickness lies beyond the domain of the validity of the thin plate model, in which case a thick plate model is required.

4.3.3. Viscous Layer Model

Once we introduce viscosity, potential theory cannot be applied, and the situation becomes significantly more complicated. Early attempts simply considered waves on a deep, viscous liquid and applied the classical solution of [71], then a more realistic two-fluid model was developed [72,73], wherein the upper layer had lower density but much higher viscosity. Various other refinements have been added to the viscous layer model but most of them share the undesirable property that the dispersion relation is not expressible in a separable form $\omega =\mathsf{\Omega }\left(\kappa \right)$ and must be solved numerically for the selected values of the dependent variable. Nevertheless, such models are of abiding importance.

To illustrate the broad features of the dispersion relation in this case, consider the Keller model for waves long compared with the viscous layer thickness,

$$\frac{g\kappa \tanh \left(\kappa H\right)}{{\omega }^2}=\frac{1-\frac{{\kappa }^{2}gh\left({\omega }^2-4i{\kappa }^2\omega \nu \right)}{{\omega }^4+16{\kappa }^4{\omega }^2{\nu }^2}}{1-\frac{{\omega }^{2}h{\rho }_{ice}}{g\rho }+\frac{{\kappa }^{2}gh\left(\frac{{\rho }_{ice}-\rho }{\rho }\right)\left({\omega }^2-4i{\kappa }^2\omega \nu \right)}{{\omega }^4+16{\kappa }^4{\omega }^2{\nu }^2}}$$

(42)

Closed form solutions for the real and imaginary parts of the wavenumber of the dominant low frequency wave mode in deep water can be obtained by setting $\left(\frac{{\rho }_{ice}-\rho }{\rho }\right)=0$ and $\kappa H=\infty$. Defining dimensionless variables for clarity of the results,

$$\tilde{\kappa }={\left(\frac{{\nu }^2}{g}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\kappa and \tilde{\omega }={\left(\frac{\nu }{g^2}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\omega$$

We obtain

$$\tilde{\kappa }={\tilde{\kappa }}_{real}+{\tilde{\kappa }}_{imag}=\left[{\tilde{\omega }}^2+2^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\tilde{\omega }}^{\raisebox{1ex}{$13$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-36{\tilde{\omega }}^8\right]+i\left[4{\tilde{\omega }}^2-2^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\tilde{\omega }}^{\raisebox{1ex}{$13$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]$$

(43)

These are plotted in Figure 18.

Somewhat surprisingly, we observe in Figure 18a that the real part of the wavenumber is almost independent of viscosity over a wide range of physically sensible values, changing imperceptibly from the free surface value. For very long periods, finite depth effects will gradually reduce the real part of the wavenumber. This warns us that we are unlikely to be successful if we try to distinguish between a viscous layer situation and a free surface by relying on the first-order echoes. For that we will need to exploit the consequences of the frequency-dependent attenuation arising from the imaginary part of the wavenumber, and plotted in Figure 18b, though these may not be evident in a resolution cell close to the ice edge, where differential attenuation has not had time (or rather, space) to modify the ice-coupled wave spectrum.

4.3.4. Mass Loading Model with Viscosity

Early attempts to model pancake ice argued that, at the early stages of development when the spatial density is low, the pancakes may act as independent masses, so the mass loading model would seem appropriate, yet experiments revealed that this model over-estimates ice thickness [74]. Subsequent comparisons with it became apparent that the viscous medium in which the pancakes formed can seldom be ignored, so a composite or hybrid model is necessary. In [75], it was demonstrated that the Keller viscous layer model gives better results for frazil and pancake ice, though of course the assumption of a homogeneous viscous layer fails to account for the mechanical interaction of pancakes embedded in the ice layer.

4.3.5. Viscoelastic Models

To overcome the deficiencies of the simpler models, attention has focused on accommodating both elasticity and viscosity in the so-called viscoelastic models, constructed so as to allow for a continuous transition between ice species by variation of the parameters. Three variants have been explored in detail. The Wang–Shen (WS) model [76,77] treats the ice as an incompressible, viscoelastic fluid in which the stress-strain relationship is modeled by a Kelvin–Voigt element. The model integrates the two mechanisms through the artifice of modifying the viscosity by addition of a term proportional to the shear modulus, yielding the effective viscosity,

$${\nu }_e=\nu +\frac{iG}{{\rho }_{ice}\omega }$$

(44)

and the implicit dispersion relation takes the form

$${\omega }^2=Q_{ve} g\kappa \tanh \kappa H$$

(45)

with

$$\begin{gathered} Q_{ve}=1+\frac{{\rho }_{ice}}{\rho }\left[\left(g^2{\kappa }^2-N^4-16{\kappa }^6{\alpha }^2{\nu }_e{}^4\right)S_{\kappa }{S}_{\alpha }-8{\kappa }^3\alpha {\nu }_e{}^2{N}^2\left(C_{\kappa }{C}_{\alpha }-1\right)\right] \\ \times {\left[g\kappa \left(4{\kappa }^3\alpha {\nu }_e{}^2{S}_{\kappa }{C}_{\alpha }+N^2{S}_{\alpha }{C}_{\kappa }-g\kappa S_{\kappa }{S}_{\alpha }\right)\right]}^{-1} \end{gathered}$$

(46)

where $S_{\kappa }=\sinh \kappa h$, $C_{\kappa }=\cosh \kappa h$, $C_{\alpha }=\sinh \alpha h$ and $N=\omega +2i{\kappa }^2{\nu }_e$.

The family of wave modes associated with the multiple roots of this dispersion relation is described in [78]. They include evanescent, gravity, shear, pressure, interfacial, and Rayleigh–Lamb waves. A further complication of such models has been proposed in [79] where the property of poroelasticity is added to accommodate ice concentration, with a view to widening the domain of the applicability of the model.

The Fox–Squires (FS) model [80] starts from thin beam theory and introduces viscosity via a complex-valued shear modulus, $G_{\nu }=G-i\omega {\rho }_{ice}\xi$, where the additional viscosity term $\xi$ is associated with the so-called ‘dash-pot’ parameter [81], leading to the dispersion relation

$$\left[\frac{G_{\nu }{h}^3}{6} \left(1+\nu \right){\kappa }^4-{\omega }^2\rho h+\rho g \right]g\kappa \tanh \left(\kappa H\right)-{\omega }^2\rho =0$$

(47)

Another variant is the Robinson–Palmer (RP) model [82]; in this case the viscosity is introduced via a frequency-dependent friction term:

$$\left[1+\frac{G{h}^3}{6} \left(1+\nu \right){\kappa }^4-{\omega }^2\rho h+\rho g-i\omega \gamma \right]-{\omega }^2\rho =0 \left(z=0\right)$$

(48)

As stated in [83,84], all three viscoelastic models can be expressed in a common form,

$$D_j\left[\left(1+Q_j\right)g\kappa \tanh \left(\kappa H\right)-{\omega }^2\right]=0$$

(49)

where the index j identifies the WS, FS or RP model.

One difficulty we encounter when implementing these models, especially WS, is the proliferation of solutions. As we are concerned with forward-propagating waves subject to frictional processes, we focus on solutions lying in the first quadrant of the complex $\kappa$ space. Many lie far from the real axis and hence correspond to rapidly decaying waves, so they are not of interest to our endeavor. It would be useful to have a means for obtaining a reasonably high-fidelity analytic expression in explicit form but one has yet to be found, however, the FS and RP models can be solved numerically with ease, whereas the WS model involves a more demanding procedure. What is particularly significant from a practical viewpoint is the fact that quite high powers of some parameters are involved, implying a corresponding sensitivity that could be exploited in measurements.

4.3.6. Discrete Floe Models

We have already mentioned the challenges of attempting to characterize the dynamics of a distribution of floes large enough to flex and otherwise involve mechanical deformations, while, at the same time, being able to collide and couple with each other and scatter waves propagating in their enveloping fluid environment. This complexity has not daunted researchers from developing models of the responses of individual floes to impinging waves, floe–floe collisions, rafting and overwash [85,86,87,88,89,90,91], and then working towards a theory of wave propagation through a floe field that takes account of multiple scattering events that are consistent with the physics at the level of the individual floes [92,93,94,95,96,97,98,99,100]. Some attempts have been made to formulate propagation through an assembly of floes as a Boltzmann equation [101,102], with the collision term now responsible for describing the hydro-mechanical dynamics of wave–floe interactions, but we understand that the jury is still out on this.

A number of experimental studies of wave propagation have been conducted to support the theoretical investigations of wave transformation during the passage through sea ice, often involving the deployment of sensors dispersed across large distances. The principal observable in most of these campaigns was attenuation, which is in accordance with the observation made earlier, that this is the more sensitive signature of changes to the complex wavenumber. One study concluded that frictional losses, not scattering, dominate attenuation mechanisms for waves whose periods exceed about 19 s [99]. This is not necessarily a positive outcome from the radar perspective because, under most circumstances, HF radar extracts more detailed information from shorter waves. More generally, unless the space–time correlation function of the surface over each resolution cell possesses a spectral representation, specifically in terms of an expansion into plane waves, it is not clear that HF radar will be able to extract detailed information about the ice field. Some measures of velocity spread, certainly, and amplitude information in the form of integrals such as significant wave height, probably, but this is not the kind of detail that we entertain as the prospect when a well-defined dispersion relation governs the dynamics.

For now, in our HF radar study, we resign ourselves to confining our attention to forms of sea ice other than fields of floes; this is the same kind of restriction employed by [103], which presents an informative selection of figures showing the frequency dependence of the resulting ice-coupled wave properties.

4.4. Model Validation and Calibration

The obvious approach to testing the fidelity of a model and thence mapping its domain of validity is to exercise it with experimental wave data where most parameters can be measured, tuning the remainder to find the best fit and then checking whether the resulting parameter values are sensible. Often this method has been adapted to compare two or more models on the same data, partially eliminating the need to solve for parameters that are common to all models. Well-known examples of this Darwinian approach include (i) comparisons of the mass loading and elastic plate models [104], (ii) the mass loading and Keller models [75], (iii) the Lamb–Stokes single-viscous layer and Keller two-viscous layer models, [105], (iv) the Keller, De Carolis and Desiderio, and De Santi and Olla viscous models [18], and (v) the WS, FS, and RP viscoelastic models [83]. Sometimes elaborate datasets are needed: the viscoelastic models are designed to handle a variety of ice types, so the first step is to establish the relation between the rheological parameters and the type of ice. An example of this procedure can be found in [106], where the WS model was fitted to wave tank measurements in which three different types of ice were prepared. The results showed an increase in effective viscosity as one moved from frazil -pancake ice, through pancake ice, to fragmented ice floes, with equivalent elasticity following the same trend.

The cost, complexity, and hazards of data acquisition campaigns in the ice zone have limited the quantity and quality of in situ measurements that can be used to validate and calibrate models, though imaginative experimental design and optimum disposition of sensors has yielded some valuable datasets [21,100,107,108,109,110]. Nevertheless, a great reliance has been placed on remotely sensed data, especially SAR. In most studies of model fidelity and calibration, researchers have taken advantage of the behavior shown in Figure 18, namely the extreme sensitivity of attenuation, associated with the imaginary component of the wavenumber, compared with wavelength stretching or compression, related to the real component. This makes excellent sense but existing technologies for remote sensing of these parameters are not entirely free of complications. For example, SAR snapshots showing a variation in wavelength across a distance of tens of kilometers may genuinely show a single harmonic wave manifesting changes in its dispersion relation because of ice changes, but it could equally show a dispersed wave packet whose leading-edge component frequencies are different from those trailing. We occasionally see this effect in HF radar echoes from the open sea when a wave packet of long waves travels through a radar resolution cell modulating the echoes from the shorter Bragg-resonant waves as it passes.

As explained in Section 3, HF radar extracts information about the full directional wave spectrum, not just the characteristics of some chosen dominant wave. This detailed information can add an extra dimension to the validation of ice-coupled wave models, and is one of the unique capabilities that HF radar brings to ice monitoring, for although SAR has some limited ability to measure the directional spectrum, it does not compare in terms of detail and accuracy. Accordingly, it is of crucial interest to understand how the different models predict the evolution of the directional spectrum of ice-coupled waves as they penetrate the ice field. This subject has certainly not been overlooked by ice researchers [80,111,112,113], and, as expected, the results are model dependent. For effective medium models, the only mechanism available for changing the directional wave spectrum in the course of propagation within the ice field is differential absorption, whereas the discrete floe models predict the redistribution of energy in the azimuth domain as a consequence of wave-floe interactions. From the radar perspective, one of the conclusions of [113] is of particular interest: “The model shows that the directional wave field in runs A and B spreads quickly and becomes isotropic or nearly isotropic in the range of wave periods T ≤ 10 s.” If this holds in the field, there will be a richer supply of waves available to support the multiple scattering and hence retrieval of the directional wave spectrum, but only as far as the short waves penetrate the ice field. This might seem to frustrate the advantage conferred by isotropy, but our modelling shows that even modest amplitudes can lead to appreciable echoes via second-order scattering when the other partner wave is strong.

5. HF Radar Signatures of Ice-Coupled Waves

5.1. Computation of Radar Signatures

With just two exceptions [34,114], where water depth was taken into account, previous modelling of HF scattering from the ocean surface has dealt with a free surface in deep water where the classic dispersion relation ${\omega }^2=g|\overrightarrow{\kappa }|$ applies, so the expressions for the various kernels in the scattering integrals are customarily found with this assumption of being ‘hard-wired’. The evaluation of the integrals can be carried out in the wavenumber domain or the frequency domain, connected by the appropriate Jacobian as determined by the dispersion relation. From the more general perspective, it is preferable to work in the wavenumber domain where the Bragg scattering mechanism retains its intuitive physical interpretation, which is the manner we have undertaken in Section 3.

The integrals for the second-order contributions must be computed along arcs in wavenumber space, as pictured in Figure 19 for the free surface case. The heavy overlaid arcs shown in black define the support of the integration.

When we apply the scattering theory to ice-coupled waves, the paths of integration shift, resulting in a modified Doppler spectrum. Figure 20 shows the revised paths for the mass loading model computed for three different mean ice thicknesses, superimposed on the free surface case. For clarity, only one quadrant of the full integration space is shown. As one might expect, the variation accelerates as ice thickness increases.

5.2. Discrimination between Ice Species

The radar Doppler spectrum depends on both the type of ice present and the parameter values for that ice species, as well as the parameters that define the interrogating radar signal. As with all inverse problems, the retrieval of the ice information we seek must contend with questions of existence, uniqueness, and stability. Existence will be an issue if no model in our library can be found to fit the data within the acceptable tolerance. Uniqueness is a serious threat because there are generally several model parameters, so it is possible for the dispersion surfaces of two different species of ice to intersect for specific values of those parameters. Stability must be dealt with by incorporating ancillary information in the retrieval procedures, typically in the form of constraints such as smoothness.

The most rudimentary radar measurement involves measuring the Doppler shifts associated with the first-order scatter. Figure 21 overlays Doppler spectra computed for two models: a free surface and a mass-loading model. For each model, the two dominant peaks, associated with advancing and retreating waves at the resonant wavenumber, have different Doppler shifts, meaning that their dispersion surfaces are separate at that wavenumber.

As evidenced by Figure 21, the high Doppler resolution of HF radar makes discrimination straightforward in the case shown, but single radar measurements of the first-order scatter can founder on the uniqueness problem. To overcome this, we can map a subset of the dispersion surface by conducting a series of measurements at different radar frequencies, so that even if one measurement lies on the intersection of two model dispersion surfaces, another may lie far enough away to enable discrimination, even allowing for a range of parameter values for each ice species. Thus, it is hard to imagine an ice monitoring radar achieving reliable classification to type without frequency agility. However, there is another pathway to discrimination, namely the second-order parts of the Doppler spectrum. The double-bounce nature of the electromagnetic scattering mechanism means that waves far apart in $\overrightarrow{\kappa }$-space are simultaneously engaged in forming the radar echoes, so the entire dispersion surface is being sampled, albeit with non-uniform weighting.

To illustrate this point, Figure 22 compares radar Doppler spectra for two models: a free surface and a mass-loaded surface. The upper figure shows spectra for a free surface, computed for three radar frequencies, while the lower figure shows the corresponding spectra for a surface covered with small ice floes modelled as mass-loading. Instead of using Doppler shift as the abscissa, we have converted to velocity. If we first look at the blue curves, corresponding to a radar frequency of 24 MHz, we see that the first-order peaks in upper and lower figures do not quite line up, reaffirming the idea that this most basic measurement can discriminate between a free surface and one with a particular ice structure that determines its Doppler frequency displacement. However, suppose that the difference was too small to measure. Now look at the second-order parts of the spectrum, extending to the edge of the plots. There are dramatic differences between upper and lower plots, revealing a powerful discrimination capability. As we go to the lower frequencies, 16 MHz and 8 MHz, the differences decrease slightly, but in the ice field it is the longer ice-coupled waves that penetrate furthest, and these respond more to the lower radar frequencies, so, even apart from radar propagation issues, there is a trade-off between diagnostic power and effective coverage.

5.3. Sensitivity to Parameter Values

Being able to map ice type synoptically and persistently may be a very useful capability but we might wish for more; we would like to be able to map the values of the ice parameters for a given species. This objective can be seen as merely a refinement of determining the ice species as discussed in Section 5.2, though with more demanding tolerances, but we can improve upon this if we use the model details to construct a measurement strategy. To illustrate how this information could help, let $B_{\alpha }\equiv \left\{{\beta }_1,{\beta }_2, \dots {\beta }_M\right\}$ denote the values of the set of $M$ parameters describing a given species of ice, $\alpha$ say, and let $C\equiv$ $\left\{{\chi }_1,{\chi }_2, \dots {\chi }_N\right\}$ denote the set of radar parameters under the control of the radar operator. Often the key radar parameter is the carrier frequency, ${\chi }_1=f$ which is restricted to some band $F\equiv \left[f_{min}, f_{max}\right]$ by both hardware and environmental conditions. We could define the ${\beta }_j$—sensitivity ${\psi }_j$ of the radar Doppler spectrum $D\left(\omega |\alpha ,B_{\alpha },{\chi }_1\right)$via the known analytical expression for the appropriate partial derivative of the identified model, together with the best available estimates for $B_{\alpha }$,

$${\psi }_j\left[D\left(\omega |\alpha ,B_{\alpha },{\chi }_1\right)\right]=\left|\frac{\partial D\left(\omega |\alpha ,B_{\alpha },{\chi }_1\right)}{\partial {\beta }_j}\right|$$

(50)

and then adopt a strategy that tunes the radar frequency to maximize the sensitivity for that model,

$${\chi }_1=f=ma{x}_{f\in F}{\psi }_j\left[D\left(\omega |\alpha ,B_{\alpha },f\right)\right]$$

(51)

Thus a priori knowledge of radar species could be used to enhance parameter measurement accuracy. This sounds rather trivial, but with HF radar, it is sometimes the case that the dimensionality of the ‘tuning space’ N is no longer small and the surface of the penalty function to be minimized has a complicated topography. For this purpose, an accelerated genetic algorithm is sometimes employed.

6. Synoptic-Scale Structure of Ice-Coupled Wave Fields

6.1. The Action Balance Equation

The mathematical models of ice-coupled waves described in the preceding sections describe the local properties of these waves, but in practical applications we must deal with extensive ice fields, where a large-scale or global representation is needed. The situation is analogous with that faced in operational oceanography, which is concerned with the evolution of ocean wave fields over large times and distances. A key tool in the development of forecasting models in that context is the action balance equation, which describes the evolution of the generalized action,

$$N\left(\overrightarrow{\kappa }\right)=\frac{S\left(\overrightarrow{\kappa }\right)}{\mathsf{\Omega }\left(\overrightarrow{\kappa }\right)}$$

(52)

due to the physical processes of wave generation by wind, nonlinear interactions that exchange energy between primary waves, and wave dissipation due to various loss mechanisms.

$$\frac{\partial N}{\partial t}+ {\nabla }_x·\left({\nabla }_{\kappa }\mathsf{\Omega }\left(\kappa \right)N\left(\kappa \right)\right)- {\nabla }_{\kappa }·\left({\nabla }_x\mathsf{\Omega }\left(\kappa \right)N\left(\kappa \right)\right)= S_{in}+ S_{nl}+ S_{dis}$$

(53)

where the source terms $S_{in}$, $S_{nl}$ and $S_{dis}$ represent the generation, nonlinear interaction, and dissipation processes.

The large-scale evolution of a field of ice-coupled waves can similarly be described by an action balance equation but, while the kinematic terms forming the left-hand side of the equation remain the same, the source terms are quite different. As in the free surface case, the solution at any given location depends on the global structure of the ice field. Some recent studies of attenuation of ice-coupled waves have embedded wave coupling models in the framework of the action balance equation [115,116], while developers continue to evaluate wave-ice parametrization models within the WAVEWATCH III forecasting model [117]. A key part of this exercise is analysis of the sensitivity of models to ice parameters [118], though few studies focus on this.

In order to specify the source terms, we need to identify the physical processes that are involved in launching ice-coupled waves, controlling their interaction with each other, and dissipating wave energy within the ice field. Here we make some observations from the radar perspective.

(i)

$S_{in}$

Although energy transfer to the ice-coupled waves from the wind above the ice surface is possible when there are substantial areas of open sea water (leads and polynyas), where sufficient fetch exists for wave growth, this type of contribution to the wave field is seldom as significant as wave energy that has propagated into the ice field from the open sea. Therefore, we might start with the model of a directional wave spectrum in the open sea in front of the ice edge and compute the wave spectrum that is transmitted through the boundary, to serve as the initial solution which can be marched deeper into the ice field, ignoring further transfer of wind energy to ice-coupled wave energy. A key part of this is the frequency-dependent transmission coefficient and frequency-dependent refraction of the directional wave spectrum. HF radar may be able to measure these.

(ii)

$S_{nl}$

Nonlinear interactions operate on the shortest timescale, and these dominate the free surface ocean wave field evolution in many circumstances. That said, triad interactions of free surface waves—the nonlinear coupling of three waves—cannot transfer energy between modes because the wave resulting from the resonant interaction of any pair of waves does not lie on the mass shell, that is, it does not satisfy the first-order wave dispersion relation. If nonlinear interactions are to be non-negligible in this scenario, they must occur at higher orders where we can simultaneously satisfy

$${\displaystyle \sum }_m{\overrightarrow{\kappa }}_m = 0$$

(54)

$$\mathsf{\Omega }\left({\displaystyle \sum }_m{\overrightarrow{\kappa }}_m\right) = {\displaystyle \sum }_m\mathsf{\Omega }\left({\overrightarrow{\kappa }}_m\right)$$

(55)

which we can interpret as conservation of momentum and conservation of energy. We note that the contribution of higher order terms decreases rapidly with increasing order except under very rapidly changing meteorological conditions.

When we consider the case of ice-coupled waves, the immediate question is whether the dispersion relation of the prevailing ice species can support resonant triad interactions. If so, then we would expect $S_{nl}$ to dominate. Otherwise, the fourth-order interactions govern energy exchange and, stretching our intuition from the free surface case, we would expect these to be weak except during sudden changes in wind stress. However, even though the triad interactions may be virtual, that is, non-resonant, HF radar has the remarkable ability to measure them because they manifest their presence as bound harmonics on the primary waves and contribute a special signature in the Doppler domain. This is exploited to great effect in the free surface situation; it may be that there is a valuable signature for HF radar here.

(iii)

$S_{dis}$

Ice-coupled waves experience attenuation when the medium is characterized by viscosity, inelastic floe collisions, flexural heating, or various other mechanisms depending on ice type. As we have noted in Section 4, the propagation loss mechanisms, associated with the imaginary component of the wave vector, emerge naturally within the effective medium treatment, primarily through the viscosity terms in young ice fields. The most important issue from the radar perspective is the rate at which the shorter waves—wavelengths shorter than 50 m, say—attenuate, because these are the most fruitful bearers of information to the radar.

6.2. HF Radar Estimation of the Source Terms

Given the parametric solutions to the dispersion relations, and exploiting the wide area coverage of HF radar to ‘follow’ the ice-coupled wave field as it penetrates the ice zone, we can envisage a unique means for studying the source terms in the action balance equation, at the same time as carrying out a diagnosis of the ice state from the local properties retrieved from the Doppler spectrum. We would expect the dissipation term $S_{dis}$ to be relatively easily retrieved, and if the radar sensitivity is sufficiently high and shorter waves present are in adequate strength, the second-order scattering terms could yield information on the nonlinear coupling term $S_{nl}$. Being able to measure the wave field in the ice over several hundred kilometers simultaneously and persistently should yield insights into the dissipation term that is inaccessible to the in situ measurement campaigns that are limited to tens of kilometers.

6.3. Ice Field Boundaries

In the preceding paragraphs we have mentioned the role of the ice boundary as the effective ‘source’ of the penetrating waves. We are concerned only with edges exposed to the open sea; shore-fast ice anchored to the shore or sea floor does not move with winds or currents, and any response to ice-coupled waves propagating through the ice field to the shore would likely be damped and highly site-specific.

Solutions describing the boundary transmittance and reflectance functions for free surface waves incident on the ice zone from the open sea have been reported by various authors [68,69,77,119,120]; some of these have been compiled and plotted in [77], which we reproduce in Figure 23.

We should regard these functions as cousins of the source terms in the action balance equation. As we will see below, a special class of HF signatures is associated with the reflectance.

6.4. HF Radar Signatures of Ice Boundary Reflectance

Ocean waves incident on the edge of the ice field from the open sea are partially reflected, analogous to electromagnetic waves encountering a medium with a different refractive index. There is a special aspect of this effect that bears on the radar signature. We have hitherto focused on the ice coupled waves within the ice field, but we should not neglect changes to the directional wave spectrum in front of the ice field which is now augmented with reflected wave energy travelling away from the ice edge. As the source term for wave generation by wind within the ice field is generally weak, and as the fetch for off-ice waves starts from zero, the reflected waves provide another, quasi-independent signature of the ice field. Specifically, HF radar can measure the full directional wave spectrum in front of the ice edge and compare it with that observed further away, where the reflected waves have been attenuated or redistributed by wave interactions. The reflectance is a function of the wave frequency and the angle of incidence as well as the ice thickness and mechanical properties, so the directional wave spectrum transformation is significant. To illustrate, Figure 24 shows the Doppler spectra for a spatial cell in front of the ice edge, computed with and without reflected waves. The frequency-dependent reflection coefficient used was that reported in [77], which was evaluated for 0.15 m ice thickness.

7. Radar Performance and Operational Limitations

The design of HF radars for high latitude operations and the formulation of optimum resource allocation strategies lie far beyond the scope of this paper, but it is important that we balance the promising theoretical predictions of radar capability in the ice monitoring role with some of the environmental, physical and instrumental factors that ultimately limit radar capability. To this end, we outline a few of these limiting factors in the following paragraphs and comment on their severity.

7.1. Limitations Arising from the Ionosphere

HF skywave radars at high latitudes must contend with more challenging propagation conditions than they enjoy at mid-latitudes. The multiscale structure of the ionosphere, the geomagnetic field geometry close to the magnetic pole and the plasma dynamics conspire to degrade measurements requiring the high phase stability of the radar signal. Signal processing techniques to treat such effects are well established [121,122,123], but they have their limitations in this complex environment [124].

In addition to these modulational effects, the small-scale plasma structures are themselves a strong source of echoes, and this applies to HFSWR systems as well as skywave radars, since it is impossible to confine all the HFSWR radiated energy to the surface wave mode. With HFSWR the echoes or clutter from the ionospheric irregularities will appear at ranges beyond some minimum value, which may be as low as 100 km, but lies more typically in the region of 150–250 km. Up to this distance, the echoes are pristine.

For these reasons, radar siting and orientation must be chosen carefully with respect to the region to be monitored and its ionospheric climatology, to minimize the deleterious effects of inclement ‘space weather’. Advances in signal processing and antenna design, together with radar deployments that respect the phenomenology, go a long way towards reducing these problems to an acceptable level but the exquisite sensitivity of HF skywave radar at lower latitudes will be near impossible to duplicate.

7.2. Limitations Related to the Resolving Power of HF Radar

Central to the feasibility of the proposed ice diagnostic technique is the issue of the Doppler resolution of the radar. If the radar cannot distinguish between the first-order echoes of two types of ice, then the primary avenue for discrimination is closed. Using conventional FFT-based spectral estimation techniques, the nominal resolution is determined by the coherent integration time (CIT), which is typically an order of ~10² s, though CITs up to 600 s have been employed. To illustrate the consequences of the spectral resolution limit, Figure 25 plots the differences between (i) the free surface dispersion relation and that for a mass loading model, and (ii) the free surface dispersion relation and that for a viscoelastic model, assuming a CIT of 100 s achieving a nominal resolution of 0.01 Hz.

The actual parameter values used for this illustration are not important here; they were chosen to highlight the fact that the longer the ice-coupled waves, the finer the resolution needed to discriminate between the models. In this picture, the two models are resolvable for short waves, up to perhaps 28 m wavelength, and the ice versus no ice situation is resolvable for waves up to 20 m for the mass loading case and 40 m for the viscoelastic case. As discussed in the preceding section, the effective Doppler resolution of skywave radar cannot be increased indefinitely because of signal distortion during ionospheric propagation.

Another point illustrated in the figure is the possibility of two dispersion curves crossing, so, at that particular wavenumber, the Doppler frequencies are identical. Of course, this can be dealt with by the simple expedient of changing the radar frequency, thereby shifting attention to waves with a different wavenumber.

7.3. Limitations Imposed by the Poor Spatial Resolution of HF Radar

In Section 3.1, we summarized the resolving power of HF radar in range, azimuth and Doppler. A spatial cell dimension of an order of 10¹ km is often achievable with skywave radar while HFSWR systems offer even greater resolution. In the grand scale of ice field formation and evolution, and away from the ice edge, it should not be uncommon for the ice within a cell of this size to have a reasonably homogeneous structure, making classification a sensible objective. That said, if distinct ice regimes are present, their joint signature will reflect all contributors, with the likelihood of degrading the inversion. Propagation considerations, HF channel availability, and antenna aperture constraints make it virtually impossible to increase radar spatial resolution significantly over present-day capabilities.

7.4. Attenuation of HF Surface Waves Propagating over Ice

It has frequently been stated by users of low power HFSWR systems that once signals leave the open sea and start to propagate over ice, the signal rapidly decays due to electrical losses in the dielectric medium. Indeed, this effect was used in [51] to map the edge of an encroaching ice field: ‘when the echo vanishes, you’ve reached the ice’. Such observations are physically correct, but we need to take account of the fact that almost all of these experiments were conducted at the upper reaches of the HF spectrum where the loss is greatest. In addition, there is a big difference between a commercial HFSWR system radiating a mean power of 30 W and a military grade system radiating 50 kW or more. Yet another factor to consider is the existence of guided wave modes when the ice cover is present. To illustrate, Figure 26 shows the electric field strength of the surface wave mode as a function of distance, computed for different ice thicknesses. It can be seen that the field is enhanced when ice is present, not diminished, over the first 20 km or so in this instance, after which it attenuates more quickly than for the free surface case. The enhancement is a function of radar frequency and can extend to much greater distances at lower frequencies; some hundreds of kilometers at 1 MHz. Therefore, as we expect ice-monitoring HF radars to spend the most time at low frequencies, the problem may not be so severe as reported.

Further, while it is undoubtedly true that measuring sea ice deep in the interior of an ice field is not feasible with many HFSWR systems, signals can certainly travel hundreds of kilometers over the sea before penetrating a distant ice field far enough to scatter from ice-coupled waves and return. HFSWR is able to detect ships reliably at ranges in excess of 300 km, occasionally up to 450 km.

7.5. Limitations Imposed by Multiple Scattering of HFSWR Signals

The notion that the signal radiated by an HFSWR propagates to a distant range cell without being distorted in the course of its passage across the intervening sea and ice field is simplistic, however, to a surprising extent, it is realized in practice. Although the perils of ionospheric modulation are avoided, multiple scattering from surface roughness and variations in ice electrical properties has the potential to distort the signal structure on both outbound and inbound paths [125]. We illustrate this possibility in Figure 27, where two signal paths are shown from target cell to receiver.

The observation that young ice fields tend to have very small surface roughness works in our favor here; multiple scattering should not contribute much in these circumstances, though quantitative studies are needed. Once the ice field matures and breaks up under wave action, floe collisions can greatly increase the surface roughness via processes of rafting and overturning. In these circumstances, the radar theory may need to be modified as it has been for the corresponding case of rough seas in ice-free waters.

7.6. Limitations Associated with the Radar Operating Frequency Band

The HF band is officially designated to extend from 3 to 30 MHz, which corresponds to first-order resonant scattering from free surface waves with periods in the range of 1.8 to 5.7 Hz. These are much shorter waves than those observed deep within the ice field. If we encroach on the MF band, operating at 1 MHz, say, the resonant waves have a period of 9.8 s, which is more promising. However, we can improve upon this. The HF radars illustrated in Figure 2 and almost all the radars pictured in Section 3.4 are monostatic or quasi-monostatic, which simplifies many aspects of radar design and minimizes cost. Bistatic radars [126], with their component sites far apart, offer a significant advantage from the ice sensing perspective: the resonant scattering mechanisms that provide the information we seek sample longer waves than those that contribute for monostatic geometries. We illustrate that in Figure 28, which maps the first-order resonant wavelength against the bistatic angle and radar frequency for an ice cover described by the mass loading model with an equivalent ice thickness of 0.5 m. Now the lowest period observable via first-order resonance reaches 15 s at a bistatic angle of 140°.

Of course, HF radar can obtain eminently usable Doppler spectra via the double scatter mechanism, sometimes even when the first-order returns are absent altogether, but the immediacy of the information in the Bragg peaks encourages us to exploit them whenever possible. A corollary to studying longer resonant wavelengths is the reduction in the associated wave frequencies and hence the need for longer coherent integration times in the observing radar.

We have mentioned the advantages of multi-frequency operation and now bistatic operation. A very informative way of visualizing the relative amounts of information that can be obtained by these different observational degrees of freedom is presented in Figure 29. Figure 29a marks the two wavevectors of the directional spectrum of the ice-coupled waves that are sampled by the first-order Bragg scatter for a single radar frequency. On this diagram, the radar is located at the 180° position. When we provide the radar to operate at a number of frequencies, we gain the extra information marked in Figure 29b. Now suppose there is another radar system interrogating the same region from a different direction, 210° here, the so-called stereoscopic configuration [126]. This adds the points shown in blue in Figure 29c. If we then allow either (or both) of these radars to receive echoes from the transmissions of the other radar, we add the points shown in violet in Figure 29d. Note how the points move to lower wavenumbers: this is the bistatic advantage in action. In Figure 29e we show the results of the joint inversion of first- and second-order echoes when a single radar is employed and signal quality is high. The lone radar is unable to resolve the left-right symmetry and hence must provide the user with both solutions. Figure 29f shows that the ambiguity can be resolved when two radars combine their measurements.

8. Opportunities with Present-Day Radars

HF radar is presently undergoing a resurgence, with a number of new systems being deployed in locations from which they can view the polar ice zones. It may happen that clients of these radars see value in developing an ice-monitoring capability. Equally, we may anticipate upgrades to existing systems that might enable them to measure ice-coupled waves. In this section, we examine some of these potential opportunities.

8.1. SuperDARN

Both polar regions are illuminated by networks of HF radars collectively known as SuperDARN, but these essentially scientific instruments are designed to measure dynamical processes in the ionosphere and thence the magnetosphere above and the atmosphere below [127,128]. Figure 30 indicate the extensive, nominal overlapping line-of-sight SuperDARN coverage zones for the two polar regions.

In terms of power and spatial resolution, these radars are, in principle, capable of measuring some ice characteristics, but the waveforms they use to study ionospheric convection, meteor trails, energy fluxes in the magnetosphere, ULF waves, tides, polar patches and many other plasma phenomena are unsuitable for monitoring sea ice. Thus, ice monitoring requires a separate monitoring mission, rather than piggybacking on the primary tasks. This has been performed in a study that investigated the use of SuperDARN radars to measure ocean surface winds and waves [129]. The results were of limited fidelity and highly variable coverage, but more advanced radar operating and signal processing techniques could go some way to improving the outcomes. A more fundamental limitation is the low radiated power, typically ~3 × 10² W mean power, compared with ~3 × 10⁵ W for a military-grade skywave radar. Yet another constraint is the restriction to radar frequencies above 8 MHz but, to balance this, a highly relevant and potentially invaluable capability has recently been demonstrated [130]: the joint operation of two SuperDARN radars in a bistatic mode, with ground echoes detectable among those from the ionosphere.

Figure 31 shows one of the radars: the Dome C radar in Antarctica.

8.2. Military Skywave Radars

Harking back to Cold War deployments, both the US and Russia have new north-looking skywave OTHR systems either in operation or under development. In all these cases, the hardware is clearly capable of conducting the measurements required for ice echo classification, so it may eventuate that this mission is eventually incorporated in the software. Canada is testing a sophisticated experimental system which may evolve into an operational system [131]. Unlike all the other radars shown in this paper, this radar employs a two-dimensional array for greater ability to suppress unwanted clutter echoes. Figure 32 shows the receiving system.

One section of the Russian Konteyner OTHR system at Kovylkino [132] is shown in Figure 33; as there are expected to be several such radars with potential for overlapping coverage over the proximate Arctic, these radars would be ideal for testing and developing ice monitoring capabilities.

8.3. Civilian HFSWR

As mentioned in Section 3.4, CODAR SeaSondes have been used quite extensively in the Arctic and the Antarctic, but so far without attempting to study ice characteristics via the dynamic signatures of ice-coupled waves. Although limited in power and resolution, the popularity of the CODAR radars [25] suggests that there may be valuable opportunities to test the ice signature theory in some future commercial deployments, noting that signal processing takes a different approach from that employed in beam-forming radars. Similar opportunities may arise with the Helzel Messtechnik WERA radars [26], though we are unaware of any present or planned deployments of these systems within the polar zones. Several other countries possess indigenous HFSWR systems and an interest in the polar regions, so we may see China, for example, exploring the application.

8.4. Military HFSWR

According to press releases, Russia intends to deploy advanced Podsulnukh HFSWR systems across its northern coastline where they could monitor the Northern Sea Route. If these systems materialize, their high power and spatial resolution indicate that they would be ideal for testing the ideas presented in this paper. The Canadian HFSWR at Cape Race ceased operation in 2007 but an upgraded system is under development [133]. (This new radar is being designed and operated by industry and academia, but most funding appears to come from the Canadian Department of Defense, so we are treating it as a military radar).

Thus, there would seem to be an adequate supply of high latitude HF radars in operation or in planning that could be used to demonstrate at least the basic concepts described in this paper, despite their suboptimum design for the ice monitoring mission. Success in a demonstration is a sine qua non for attracting the level of interest and investment needed to develop an operational capability.

9. Conclusions

The possibility of measuring sea ice type and parameters with HF radar via the dynamics induced by ice-coupled waves originating in the open sea is intriguing and could have substantial practical benefits to activities in the polar ice zones. Our investigation has tried to approach this question with an open mind; if any undue optimism or pessimism is detected by the reader, it is a symptom of many years of experience with the extraordinary capabilities of HF radars in other missions.

It was observed in [83] that “Simultaneous measurements of wave period, attenuation rate and wavelength are required to perform calibration, which has never been achieved and may not be feasible with contemporary measurement techniques.” To the extent to which HF radar can be regarded as a contemporary measurement technique, we beg to differ, as these are precisely the observables accessible to HF radar measurement (at least for a sophisticated radar system). Under favorable conditions, for waves whose period does not exceed 10 s, the wave frequency may be measured to within 0.01 Hz, the wavelength to better than 1% (depending on the spatial resolution desired), and the attenuation rate to perhaps 0.2 dB/km. Longer waves are still observable, but radar measurement accuracy will decrease. There are other caveats that might be applied to these numbers, but they relate more to the design and siting of the radar and the prevailing space weather (for skywave radars) than to the physics of the radar measurements.

The present study is exploratory in nature, and many important questions remain to be addressed, but the broad conclusions that we have drawn may be summarized as follows:

• The persistent mapping of the ice edge on a synoptic scale with a resolution of order 10 km and a resampling interval of order tens of minutes is eminently achievable with HF skywave radar;
• HFSWR can achieve even finer spatial resolution but at reduced ranges, out to hundreds of kilometers instead of thousands;
• Unlike space-borne SAR, which relies on a snapshot of ice waves to estimate the dispersion, HF radar measures both spatial dimension AND phase speed of the waves, yielding far higher discrimination. Moreover, HF radar can achieve its analysis even in a confused sea when many waves are superimposed, where SAR might fail to identify a candidate wave for measurement;
• The dispersion relations of different types of ice support a crude classification capability using the first-order radar echoes from a radar cell under test, though not necessarily for all ice species. This result opens the door to many low-cost radars interested in participating in the ice monitoring mission, noting that radar design impacts on the retrieval;
• For a given ice species, a similar approach can yield estimates of dominant parameters, though not necessarily for all ice species;
• When the second-order Doppler spectra are included in the echo analysis, much greater discrimination between ice species and between parameter values ensues, but the form of the directional wave spectrum deep within the ice field has a strong bearing on this;
• The evolution of the ice-coupled wave spectrum as it penetrates deeper into the ice field enables the radar to take advantage of the frequency dependence of the attenuation due to viscosity, thereby enhancing ice classification, and also contributing to our understanding of the source terms in the action balance equation;
• Additional radar signatures of the ice characteristics may be found in front of the ice field, as a result of waves reflecting from the ice edge;
• Nature has decreed that the waves of greatest relevance deep within the ice field respond preferentially to low radar frequencies. This may limit the ability of some HF radars to make meaningful measurements;
• For HFSWR systems, the attenuation of the radar signal as it propagates across the ice field will be rapid once a threshold distance is passed; this distance may be just a few tens of kilometers at upper HF frequencies but can be many hundreds of kilometers at the low end of the band;
• HF skywave radars avoid this and some other problems but must deal with the high latitude ionosphere with its many irregularities and turbulent motions that can degrade signal purity and hence obscure the features needed for ice echo analysis. The efficacy of advanced signal processing techniques that work exceptionally well at lower latitudes has yet to be established;
• The morphology of multi-year ice is complex and its dynamics unlikely to follow the predictions of a simple model, but a radar signature will nevertheless appear. However, practical concern focuses on rapidly evolving first-year ice fields.

We have identified a number of areas for priority attention in future research, especially those related to (i) using existing in situ data for quantitative modelling of radar performance, (ii) inverse methods for information retrieval, (iii) techniques that optimize radar resources for the ice monitoring mission, (iv) the electromagnetics of radio wave propagation over heterogeneous ice fields, (v) techniques for accurate skywave radar coordinate registration in the Arctic, and (vi) enhancing signal decontamination procedures for the high latitude ionosphere.