Monitoring Sea Ice Thickness Using GNSS-Interferometric Reflectometry

This letter presents the analysis of frozen sea surface properties using low-cost and low-complexity terrestrial global navigation satellite system (GNSS) receivers. Monitoring sea ice thickness and the mean sea level (MSL) of the frozen sea are performed using the interference frequency obtained by the GNSS interference pattern (IP) technique. The height variations between the GNSS antenna and the sea surface evaluated using the IP of the direct and reflected carrier-to-noise density ratio (C/N0) are used to find the corresponding MSL. The GNSS-reflectometry (GNSS-R) derived MSL for open sea conditions agreed well with the mareograph data with a root-mean-squared error (RMSE) of 2.72 cm with an R-squared value of 0.9644. For frozen sea, a notable difference was observed between the measured MSL and ground-truth MSL values. This difference was caused by the combined thickness of snow and ice above the frozen sea surface, also known as the total freeboard. Assuming the conditions for hydrostatic equilibrium is satisfied, total freeboard was converted to ice thickness. The ice thickness values agreed well with the published ice charts by the Finnish Meteorological Institute (FMI). The main uncertainty in the extracted ice thickness was due to the thick snow accumulation and unknown snow properties.


I. INTRODUCTION
G LOBAL navigation satellite system-reflectometry (GNSS-R) is a passive remote sensing technique that uses reflected GNSS signals from surface of the Earth to extract various surface properties. GNSS-R has been under constant development for more than two decades and has expanded its capability in several applications, such as measurements of sea level, snow depth, snow water equivalent, wind speed, ice thickness, soil moisture, hurricane detection, vegetation characterization, and so on. GNSS-R has evolved over the years, and several different techniques for analysis have been developed, such as delay analysis, delay-Doppler analysis, carrier phase analysis, and signal-to-noise ratio (SNR) or carrier-to-noise density ratio (C/N0) analysis. Commercial GNSS receivers that provide raw RF data are expensive, bulky, and have high data payload as compared with the receivers that output only SNR or C/N0 as the measure of signal strength. Therefore, in this study, the C/N0 method will be used to analyze the direct and reflected GNSS signals. The main objective of using this method is to reduce the computational complexity, costs, power consumption, and the weight of the measurement setup. A GNSS-interferometric reflectometry (GNSS-IR) method is used to study the interference pattern (IP) created by the combination of direct and reflected GNSS signals. Direct signal is right-hand circular polarized (RHCP), whereas the polarization of the reflected signal depends on the properties of the scattering surface. The received RHCP IP possesses the information regarding the reflecting surface. The frequency of IP is the function of the vertical height of the antenna phase center from the horizontal smooth reflecting surface, the amplitude of the IP gives the information of dielectric properties, and roughness of the reflecting surface and the change in phase of IP can be attributed to change in electromagnetic properties or geometry of the reflecting medium. The IP has been studied to measure snow depth, sea ice detection, sea ice thickness, sea state, sea level, and so on [1], [2], [3], [4], [5].
This letter presents results from long-term sea surface sensing campaign using the low-cost, low-complexity GNSS-R system. The letter focuses on the possibility of monitoring sea level and sea ice thickness simultaneously using GNSS-IR antenna height variations derived using multiple GNSS satellites. Measurement scenario and setup are described in Section II. The principle of GNSS-IR and IP analysis for sea-level monitoring is given in Section III. Extraction of total freeboard and sea ice thickness from GNSS-IR data are detailed in Section IV. The letter is concluded in Section V.

II. MEASUREMENT SETUP AND SCENARIO
A dual circular polarized antenna (DCPA) was used to receive, direct plus reflected RHCP and reflected left-hand circular polarized (LHCP) signals, simultaneously [6]. The received LHCP signal contains reflected component of the signal from the snow or sea-ice surfaces as well as from the ice-water boundary and is beyond the scope of this study. Therefore, this study will focus only on the RHCP interferometric (direct plus reflected RHCP) received signal captured by the right-hand port of the DCPA antenna. A lowcost and lightweight U-blox M8T universal serial bus (USB) GNSS receiver was used to record the GNSS signals. The This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ parameters, such as C/N0, satellite elevation and azimuth angles, and time (UTC), generated in the National Maritime Electronics Association (NMEA) format are used for analysis. The receivers were connected to a Raspberry Pi computer. The LTE USB modem was used to enable remote access to the measurement system and GNSS data. The GNSS receivers, computer, modem, and the GNSS antenna were assembled inside a commercially available, regular use 10-L bucket. The bucket acted as a radome, protecting the equipment from harsh environmental conditions, also significantly reducing the time, cost, and complexity of designing and fabricating a dedicated radome. Performed antenna pattern measurements showed that the bucket reduced the gain of the DCPA by 0.1 dB in all elevation and azimuth directions.
The GNSS-R measurement system was installed at the Port of Oulu, facing the Gulf of Bothnia in the Baltic Sea. The measurement system was mounted on a wooden pole at the height of ≈3 m from the ground and ≈6 m from the sea surface (it should be noted that the height of the antenna above sea varies with the sea-level variations). The GNSS antenna was mounted horizontally to be able to receive the direct and the reflected signals with equal gains. The antenna boresight direction (θ ant = 0 • ) was pointing ≈335 • in the azimuth direction with respect to the North. The measurement setup and the frozen Baltic Sea are shown in Fig. 1(a). The signals were only recorded for satellites passing from ±30 • in the azimuth and from 2 • to 15 • in the elevation with respect to the antenna boresight direction to ensure good quality reception of RHCP interference signal. The satellite footprint that contributes significantly to the signal reflecting in the specular direction toward the receiving antenna is defined by the first Fresnel zone (FFZ). The FFZ for the considered satellite elevation and azimuth angles is shown in Fig. 1(b) [7].
The GNSS-R measurement began during the peak of the Finnish winter season, in February 2021. The average temperature for January-February 2021 was −15 • C. The ground-truth data, such as sea level, snow depth on land, temperature, and ice thickness (ice charts), were downloaded from the Finnish Meteorological Institute's (FMI) website [8]. The Baltic Sea was completely frozen with the mean ice thickness in the range 25-50 cm according to the ice chart published by the FMI at the time of installation. The frozen sea surface conditions were also manually investigated during the spring season. An ice drill was used to manually drill holes through snow and ice to measure their respective thicknesses. The manual measurements were performed on a single day, at five locations, separated by at least 50 m and ≈1 km from the GNSS-R measurement site. The drill area is shown in Fig. 1(b) (red polygon), and the measured snow, ice thicknesses, and observed sea state are highlighted in Table I. From Table I, it is clear that the snow and ice thicknesses were very localized, and also, multiple incidents of flooding between snow and ice boundary were observed. The sea surface was also visually inspected later in the spring, and the sea ice was observed to be free of snow. The snow depth data on top of sea surface were not available; therefore, the available data for snow depth on land are scaled to match the ground-truth conditions during the drill measurement.

III. MONITORING SEA LEVEL USING GNSS-IR
The RHCP received signal is the combination of direct RHCP and reflected or scattered RHCP signals from the horizontal surface. The reflected signal travels an excess path length (EPL) relative to the direct signal, which introduces a phase difference between the signals. The continuous movement of satellite in the elevation plane causes the phase difference between the direct and the reflected signals to change continuously, resulting in a pattern of constructive and destructive interference. Hence, this phenomenon creates an IP with near constant frequency at lower elevation angles. The frequency of the IP can be related to the geometry of the signal reception, which is the function of distance of the antenna from the reflecting surface and the wavelength of the carrier signal. The geometry-based EPL for a horizontal reflector is given as in [9] EPL = 2hsinθ el (1) where h is the height of the antenna above reflecting surface, and θ el is the satellite elevation angle. The phase variation due to the EPL and the satellite motion in elevation plane is where λ is the carrier wavelength. The frequency of IP is derived by differentiating (2) with respect to sine of elevation angle [10] f IP = 1 2π It can be seen from (3) that f IP is directly proportional to the antenna height. Therefore, higher frequency of IP corresponds to larger distance between the reflecting surface and antenna phase center. A spectrum analysis of the received RHCP IP signal is performed to extract the dominant frequency, which is the function of the height of the antenna above the reflecting surface. Lomb-Scargle periodogram (LSP) will be performed as the function of sine of elevation angle to extract the dominant frequency from RHCP IP. LSP performs fast Fourier transform (FFT) for non-uniformly spaced samples. Thus, height of the antenna above the reflecting surface can be directly calculated interchanging the terms in (3) with f IP extracted using the LSP. The retrieved antenna heights using multiple satellites are then arranged according to the mean time of the measurement. Moving average of eight antenna height samples was performed to remove the noisy trend in the long time series of GNSS-IR antenna height data.
GNSS-IR frequency variations are used to keep track of the antenna height above sea surface over time. The higher the IP frequency, the larger the distance between the sea surface and the height of the GNSS antenna, which also corresponds to lower sea level and vice versa. The antenna heights are converted to mean sea level (MSL) by finding the correct reference of antenna height to the corresponding sea level [5]. The most accurate reference for MSL was achieved when the sea was free from snow and ice. The GNSS-IR antenna heights where h MSLR is the antenna height corresponding the reference MSL. h MSLR is estimated by linear fitting the GNSS-IR antenna height variations and mareograph MSL data, as shown in Fig. 2 (right) [5]. The linear fit is performed from day 72 onward, i.e., for open sea conditions marked by arrow in Fig. 2  and FMI sea-level data follow similar trend throughout the three month period. However, a notable difference is observed between the measured MSL and the FMI data for the frozen period. The snow and ice accumulation on top of sea surface has to be taken into account for the frozen period for further analysis.

A. Total Freeboard Estimation Using GNSS-IR
An illustration of the reflection of GNSS signals from frozen sea surface is shown in Fig. 3. Specular reflection of GNSS signals occurs at air-snow or snow-ice boundary, depending on Fig. 3. Reflection of GPS signal from frozen sea surface. Snow-ice freeboard also known as total freeboard is the total combined thickness of snow and ice above the water surface.
the density and thickness of snow on top of sea ice. Therefore, snow and ice accumulated on top of the sea surface causes the GNSS-IR retrieved MSL to deviate from the actual MSL. This deviation is proportional to the combined thickness of ice and snow above the sea water surface, which is known as snow-ice freeboard or total freeboard. The difference between the GNSS-IR measured MSL and the FMI MSL is used to estimate the total freeboard height h fb , using the relation where SWL is the sea water level given by FMI mareograph data, whereas SSL is the sea surface level equivalent to MSL GNSS R from (4). In the absence of snow and ice, SWL ≈ SSL and h fb = 0, which signifies open sea conditions. The total freeboard for the measurement period of 97 days is shown in Fig. 4. Moving average filtering with a window length of ten days is performed to reduce time synchronization mismatch between GNSS and FMI data. Snow depth and 24-h average air temperature based on FMI data are also plotted in Fig. 4. Here, snow depth on land is rescaled to top of sea ice according to the ground-truth conditions (refer to Section II and Table I), and assuming that the GNSS signals are reflected from the air-snow boundary. According to Fig. 4, the total freeboard during first 20 days decreases rapidly but does not show any correlation with temperature. The sharp increase in snow depth is suspected to trigger this decline in a total freeboard value. Heavy snow loading on top of sea ice can push the ice under the sea surface that allows sea water to flood the boundary between snow and ice, resulting in a decrease in the total freeboard [11]. The manual ice drill measurements on day 39 also showed flooding between snow and ice boundary on three out of five drilled locations, as shown in Table I. The snow conditions as well as total freeboard remained almost constant between days 20 and 40. The snow melted quickly as spring season arrived, and the sea ice after day 50 was free of snow. The total freeboard has zero mean value after day 72, which signifies open sea conditions.

B. Estimating Sea Ice Thickness Using Total Freeboard
Measuring ice thickness can be more challenging as compared with monitoring sea ice extent. Ice thickness and ice Fig. 4. Total freeboard as the difference between the GNSS-IR sea level and the FMI sea level (black). Range of snow depths measured during manual drill measurements is shown as error bar. Snow depth data on land provided by FMI rescaled to approximately match the ground-truth conditions at sea surface (yellow). The 24-h average air temperature for the measurement period is plotted below (bar graph). structure are very localized in nature; therefore, most measurement techniques provide ice thickness data as an average over large areas. In this study, the specular reflected GNSS signals from the sea surface are considered, and the estimated ice thickness is calculated as an average over the area covered by the FFZ, as shown in Fig. 1(b). The total freeboard is converted to ice thickness, assuming snow and ice on the sea surface satisfy the conditions for hydrostatic equilibrium, given by the relation [12] where ρ w , ρ i , and ρ s are the density of sea water, sea ice, and snow, respectively, whereas h fb and h s signify total freeboard height and snow depth, respectively. The sea ice thickness is the function of densities of snow, sea water, and ice. Due to unavailability of the local sea water density values, sea ice thickness is plotted based on range of known density values found from the literature [13], [14], [15]. In Fig. 5(a), ice thickness is plotted as a function of varying sea ice density, and constant snow and water density values, ρ i = [850, 870, 890] kg m −3 , ρ s = 400 kg m −3 , and ρ w = 1020 kg m −3 , respectively, whereas in Fig. 5(b), ice thickness is plotted by varying the snow density and keeping water and ice densities constant, ρ s = [300, 400, 500] kg m −3 as ρ w = 1020 kg m −3 and ρ i = 870 kg m −3 , respectively. Ice thickness in Fig. 5 shows the upper and lower bounds of the possible ice thickness values for the chosen sea ice and snow densities. In Fig. 5, the calculated ice thickness for the first 20 days shows large and unrealistic values. This could be due to unavailability of the accurate snow data for freeboard estimation, the different dynamics of the snow accumulation over land versus over the sea, and uncertainty of snow density on the sea surface. However, with the stable snow depth after day 20, the accuracy of GNSS-IR measured ice thickness compared with the FMI ice charts gets better. The ice thickness estimates during this period (from day 30 to day 50) are lower than the average but within the range given by the FMI ice charts. The most accurate estimation of sea ice thickness is observed for the period when the sea ice is free of snow (i.e., after day 50 in Fig. 5).

V. CONCLUSION
GNSS-IR was used to monitor sea-level variations through the winter 2020 and spring of 2021. The MSL derived from GNSS measurements agreed well with the mareograph data with RMSE 2.72 cm for open sea conditions. However, a notable difference witnessed between the GNSS-IR measured MSL and ground-truth data for the frozen period was found to be caused by the formation of sea ice and accumulation of snow on top of the sea ice, also known as total freeboard. Total freeboard was calculated by taking the difference between the GNSS derived MSL and MSL data from the FMI. A sharp decrease in total freeboard (unrelated to temperature) with the increase in snow depth in the same period agreed with the flooding observed between snow and ice during manual drill investigation. Total freeboard was converted to the corresponding ice thickness using the hydrostatic balance equation. The calculated ice thickness values agreed well with the published ice charts by the FMI. The main uncertainties in the extracted ice thickness were due to the unknown snow properties that needs more detailed study in the future. This study successfully showed that GNSS-IR measurements captured the sea level trends throughout the winter and springs season enabling good estimates of sea ice thickness and MSL and possibly detect event, such as flooding between snow and ice boundary. Monitoring the sea from the beginning of the freezing season will enable better estimation of ice and snow thickness on top of sea surface and help reduce the uncertainty observed due to unknown snow conditions. It will also enable the low-cost GNSS-R system to be an independent sea surface monitoring system.