Toward the Estimation of Oil Slick Thickness on Newly Formed Sea Ice Using C-Band Radar Backscatter

Oil thickness in oil spills involving sea ice is a key parameter required for an effective oil spill response; however, quantifying it from radar backscatter data remains a difficult task. We investigated a possible solution for estimating oil slick thickness by using electromagnetic (EM) forward and inverse scattering models of oil-covered newly formed sea ice (NI). Our forward model employs a first-order approximation of a multilayered small perturbation method (SPM) to predict two copolarization C-band radar backscatters of NI covered by an oil slick with thicknesses ranging from 0 to 7 mm. The results showed that the backscatter decreases as slick thickness increases, which we attributed to signal attenuation within the saline-oil layer. Our inverse model relies on the particle swarm optimization (PSO) algorithm to determine the slick thickness on NI using synthetic backscatter data, and it requires the input of several important ice and oil physical parameters (thickness, dielectrics, and roughness). Moreover, the estimated slick thickness was validated using scatterometer data from an oil-on-ice experiment at the University of Manitoba’s Sea-ice Environmental Research Facility (SERF). With synthetic data, the 5 mm oil slick thickness was overestimated by 25%, while with experimental data, it was overestimated by 8%. Overall, our findings have laid the groundwork for future inversion studies to identify the thickest oil spill zone from current and future C-band radar satellites for immediate response.

to such hazards in this remote and difficult to access region.Satellite-based microwave remote sensing (e.g., synthetic aperture radar, SAR) is a well-established tool that shows strong potential for the task [5], [6]; however, the complex physical interaction between the sea ice and oil makes the microwave interpretation of oil spills extremely difficult.This is where electromagnetic (EM) modeling process comes into play and it has two parts [7].The first part is forward modeling, where the material parameters are known and we seek the expected EM fields, and this allows us to understand how the physics of oil-ice interactions relate to the microwave scattering behavior.The second part is inverse modeling, where we take EM measurements and attempt to infer the material properties that caused those measurements.
Inversion could allow us to estimate oil slick thickness that may be distributed beneath, within, and on top of the ice, but it becomes particularly difficult because different combinations of physical properties of layered media (such as roughness and dielectrics) can cause similar backscatter signatures.This non-uniqueness is especially bad when considering scatterometer measurements, which can have limited scanning geometry and a single transmit/receive position (i.e., monostatic radar).
The ultimate goal of this manuscript is to estimate oil slick thickness on experimental oil-covered newly formed ice (NI) using ground truth C-band scatterometer data, which is particularly relevant for calibrating contemporary C-band radar satellites, including the Canadian RADARSAT Constellation Mission (RCM).To achieve this, we seek to answer two research questions as follows: 1) What is the expected C-band radar backscatter response when simulated oil slick thickness on NI increases?2) How feasible is it to estimate the slick thickness on NI using the experimentally measured C-band radar backscatter?The first question is addressed by employing a multilayered-SPM forward solver to simulate the total backscatter response of oil-covered nilas ice, while the second question is addressed by using an inverse solver (that calls the forward model) based on particle swarm optimization (PSO) algorithm.Synthetic backscatter data will be generated to evaluate the performance of the proposed PSO algorithm.As a validation step, we use scatterometer data from an oil-on-ice experiment at the University of Manitoba's Sea-ice Environmental Research Facility (SERF).
The remainder of this manuscript is organized as follows.Section II introduces key physical and radar aspects of oilcovered NI.Section III provides an overview of our EM modeling strategy used to estimate slick thickness on the ice surface.Section IV presents the results of the forward and inverse scattering models of oil-covered NI.Section V discusses the model sensitivity to changes in input oil parameters, along with the oil thickness retrieval.Finally, Section VI concludes with a summary and recommendations.

II. THEORETICAL BACKGROUND
To accurately determine the thickness of an oil slick on NI, it is important to first understand how microwave energy interacts with sea ice and oil as a function of their geophysical and thermodynamic properties.This section provides concise background information on oil-covered NI, including their dielectric mixture model, scattering measurement/model, and slick thickness retrieval.

A. Oil-Covered Newly Formed Sea Ice (NI)
NI is a complex multiphase material consisting mostly of ice crystals and brine inclusions.This type of sea ice is prevalent during freeze-up seasons and exists in various formations (such as frazil, grease ice, nilas, and pancake ice [8]), and its surfaces can be bare, covered with frost flowers, or covered with snow [9].Under relatively calm ocean conditions, NI usually grows into two distinct ice types: dark nilas (<5 cm thick) and light nilas (5-10 cm thick) [10].The topmost layer of these ice types is characterized by a highly saline frazil [11], making the NI profile a stratified material [12].
In the event of an oil spill on Arctic NI-covered waters, weathering processes are expected to occur at a slower rate or remain unchanged due to subzero ambient temperatures [13].This causes the spilled oil to spread and drift slowly over the ice terrain, creating a much thicker slick layer suitable for effective response methods such as in situ burning and chemical dispersants.

B. Dielectric Mixture Model of Oil-Covered NI
Using dielectric mixture models, the dielectric properties corresponding to each layer of the NI profile could be derived from the ice geophysical properties.The refractive dielectric mixture model has consistently shown agreement with direct dielectric measurements among various models used to calculate sea ice dielectrics [14].In the case of oil contamination in sea ice, the presence of oil introduces a new component within the ice structure, necessitating the inclusion of oil in the dielectric mixture model [15], [16].If the spilled oil forms a slick cover solely on the NI, it creates a heterogeneous mixture of oil and brine due to the high salinity of the topmost NI layer [17].The precise mixing process remains uncertain, as the hydrophobic nature of oil suggests that it will remain separate from the liquid brine, which is primarily water [18].As such, it is appropriate to employ two distinct constituent materials in the refractive mixture model, a technique also applicable to mixtures in the liquid phase [19].

C. Scattering Measurement of Oil-Covered NI
Sea ice backscatter refers to the signal that is returned when microwave radiation interacts with the surface and medium of the ice.Its strength is influenced by three key factors [20]: the ice surface roughness, measured by the root-mean-square (rms) height and correlation length; the ice complex dielectric constant (CDC), which depends on its geophysical properties; and the microwave system parameters, such as frequency, incidence angle, and polarization.Meanwhile, the penetration depth of the EM wave is mostly governed by the ice dielectric constant and radar frequency [21].To quantify the amount of backscatter per unit area of an ice scene, we utilize a metric called the Normalized Radar Cross Section (NRCS) [22].Calculating the ratio of this metric in VV (vertical transmit-vertical receive) to HH (horizontal transmit-horizontal receive) polarization will reduce the surface roughness effect on the backscatter response [23].
NRCS data can be obtained through direct observation using remote microwave sensors, such as SAR and scatterometer.For oil-contaminated scenes, the measured scattering responses are very low and close to the sensor's noise floor, known as noise-equivalent sigma zero (NESZ) [24].This leads to a low signal-to-noise ratio (SNR) of the NRCS values [25].According to Minchew et al. [24], at least 6 dB above the NESZ level will enhance the SNR of the measured NRCS for oil-on-open water.Similar studies have proposed thresholds ranging from 7 to 9 dB [26], and even up to 10 dB [27].However, in case of oil-on-ice, there is no detailed SNR analysis due to the lack of a major oil spill in ice-covered water regions.As a result, Johannson et al. [28] utilized two paired SAR images of oil-contaminated open water and oilfree NI, applying 2 dB threshold for the SNR analysis, while acknowledging the influence of noise on signal below the 6 dB threshold.We intend to use the same approach for our SNR analysis.

D. Scattering Model of Oil-Covered NI
Another way of obtaining NRCS data is via simulation, achieved by employing forward scattering models, which are classified into surface and volume components [29], [30].Because of the high salinity of the topmost layer of NI [12], [31], radar signals usually interact with the ice surface, making surface scattering the dominant scattering mechanism.
Several mathematical models have been developed to predict ice-scattering mechanisms dominated by surface scattering.These include the small perturbation method (SPM), Kirchhoff approximation, integral equation method, and computational EM techniques [29], [32], [33], [34].Of these, SPM is highly suitable for predicting the scattering of NI because it is computationally inexpensive and can simulate radar scattering of a slightly rough surface (that could be submerged in the medium), such that the vertical variation (rms height, h) and horizontal variation (correlation length, l) are both smaller than the incident wavelength (λ ).This means that SPM is applicable under specific roughness conditions, where h < 0.05λ , l < 0.5λ , and h/l < 0.3 [35].Even in situations where oil covers the ice and smoothens the surface, the SPM conditions still hold true.
Given the stratified nature of nilas ice, a multilayered SPM model will be useful in modeling the total scattering responses.Komarov et al. [36] developed a first-order approximation of the multilayered-SPM and used it to simulate VV and HH polarization backscatter of snow-covered sea ice systems when the dominant scattering occurs at the air/snow, snow/ice, and ice/seawater interfaces [36], [37], [38].Isleifson et al. [15] successfully applied this model to various oil-contaminated sea ice systems, with oil distributed beneath the ice.In our case, we implement this multilayered-SPM in a scenario where oil is on top of the ice.

E. Retrieval of Slick Thickness on Oil-Covered NI
Turning to how oil slick thickness on NI is retrieved, inverse models rely on optimization algorithms to minimize the objective function between the modeled and measured NRCS data [39].The process involves repeated calls to forward models during the iterative minimization.Prior to the actual inversion process, synthetic NRCS is often generated to test the effectiveness of such optimizers.Retrieving oil thickness on sea ice is extremely difficult because no operational radar for an oil spill incident in the Arctic is available.This limitation can be addressed by conducting experimental studies to collect in situ radar scattering of oil-contaminated sea ice.To date, only two studies by Firoozy et al. [40] and [41], have utilized a differential evolution optimizer to retrieve the oil layer thickness entrapped beneath and within ice.In contrast, our study employs a PSO algorithm to retrieve the oil layer thickness on top of the ice.Experimental studies in this research domain are scarce, and attempts to extract slick thickness when spilled oil spreads over sea ice scenes have not yet been fully explored, motivating our investigation.

A. Input Datasets
Table I provides the input physical datasets used to drive the EM modeling strategy and they were collected along with C-band scatterometer data during the SERF 2020 oil-on-sea ice experiment [17].Although the details of this experiment have been discussed elsewhere in [17], [42], and [43], it is worth noting that the sea ice was grown in a cylindrical tub of 4-m diameter and 1-m depth, under two phases of oil spill scenarios.In phase-1, oil was injected beneath existing ice, and in phase-2, the tub was filled with artificial seawater contaminated with 6.75 cm 3 of Tundra Crude (a crude oil type provided by Tundra Oil & Gas in Canada).This article focuses on the phase-2 experiment.It is important to note that we did not actively control the thickness of the oil slick.Instead, the oil thickness naturally formed a skim layer over the established ice pack, with an average thickness of 5-mm measured from core samples.
During the experiment, we measured the in situ physical characteristics of ice and crude oil, which is required for calculating the CDCs of the spilled oil, oil-contaminated ice, and uncontaminated ice.These dielectric values, along with other information on ice thickness, oil thickness, and roughness of the oil/ice (or air/ice in the oil-free case) and ice/water interfaces, are fed into our EM forward solver.We vary a number of inputs to achieve the purposes of our forward solver (see Section IV).
Throughout the experiment, our C-band polarimetric microwave scatterometer was mounted on a 5.3 m scaffolding tower and scanned over the tub at a 24.5 • incidence angle with respect to the nadir without interference from the edges of the tub.The measured backscatter values were VV = −39 dB, HH = −42 dB, and VH = −46 dB [17].These values remained relatively constant above the NESZ level of −45 and −50 dB for the co-and cross-polarization signals, respectively.As a result, we expect a poor SNR for our measured NRCS that could affect the retrieval parameters, including the oil thickness and dielectric property [44].It is worth noting that we used only the measured NRCS in VV and HH polarizations because our NRCS simulation model is limited to copolarization channels (see Section IV).
Our scatterometer has a very low NESZ compared to the current satellite SAR (e.g., RCM ScanSAR low noise beam mode [45]), which compelled us to use the 2 dB threshold from Johannson et al. [28] as a "mild threshold" above the NESZ for both modeled and measured NRCS data, ensuring their suitability for subsequent analysis.Despite this, we acknowledge the 6 dB threshold recommended by Minchew et al. [24] as a "harsh threshold" to avoid any possible noise contamination from the receiver noise background [24].Details on our scatterometer system and calibration methods are provided in the literature [46].

B. Electromagnetic (EM) Modeling Strategy
Fig. 1 provides an overview of the EM modeling strategy used to determine the thickness of the oil slick on oil-covered NI (i.e., the target).The procedure begins with the inputs of parameters (see Table I) into our multilayered-SPM forward solver, which simulates C-band radar with VV and HH polarizations and 20 • -60 • incidence angles.To serve the purposes of our forward solver, we vary a number of input parameters.In Fig. 1(a), we adjusted only the oil thickness from 1 to 7 mm to predict the expected NRCS of the target.In Fig. 1(b), we modified certain parameters as unknowns (e.g., oil thickness, oil dielectrics, and ice dielectrics) and fixed certain parameters as a-prior information (e.g., ice thickness, roughness, and water dielectrics) for each forward solver that our PSO-based inverse solver called.During the inversion process, the output is fed into the PSO algorithm which attempts to find the parameter of interest (oil layer thickness) by modifying the unknowns until the modeled NRCS matches the experimental NRCS.
Using synthetic NRCS (σ o(syn) in dB units) data serves as a preliminary step to test the applicability of the PSO and to fine-tune its hyperparameters.This better prepares our inverse solver for validation with ground truth scatterometer NRCS measurements of experimental oil spills on artificially developed saline ice.During the testing of the inversion process with σ o(syn) , we added 1% Gaussian white noise (G n ) to the where the subscripts p and q denote the linearly transmitted and received polarizations, which are orthogonal to each other as vertical (V) or horizontal (H) components, respectively.The notation, r represents the independent normally distributed random numbers at VV and HH polarizations.Note that the noise was added to avoid overly optimistic estimates of inversion performance on synthetic data, and the exact magnitude of that noise was a choice made by matching the noise until the synthetic data resembled the experimental data.
1) Forward Solver: This study utilizes a forward solver based on the first-order approximation of multilayered SPM to simulate the backscatter of an oil-covered NI profile.This forward solver was developed by Komarov et al. [36] and has been widely used for applications that require robust prediction of multilayered media, such as snow-covered sea ice [36], [37], [38] and oil-contaminated sea ice [15].Full details can be found in [36].
In the multilayered-SPM framework, the idealized ice profile is a multilayer medium between the upper air and lower water half-spaces.We set up this framework with two different ice profiles, as shown in Fig. 2. Each layer is characterized by specific parameters, such as thickness ( d), CDC (ε), and roughness statistics (h and l).The first case pertains to oilfree NI, which necessitates the calculation of NRCS from two rough interfaces (air/ice and ice/water).This model case Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.I.
provides a baseline for comparison with the oil-contaminated ice.In contrast, the second case involves an oil-covered NI, requiring the calculation of NRCS from two rough interfaces (oil/ice and ice/water).In this case, we varied the oil layer thickness on NI, in the range of 1-7 mm.The oil-covered NI has two significant aspects.First, we designate the slick oil on NI as "saline-oil" due to its contamination with brine that migrated from the highly saline topmost layer of NI.Second, we refer to the frazil layer as "oiled-frazil" because it has been contaminated by the spilled oil.
Table I shows the input parameter values used in our multilayered-SPM forward solver.The pure sea ice CDC, brine CDC, and brine volume were calculated using formulas developed by Mätzler and Wegmüller [22], Cox and Weeks [47], and Stogryn and Desargent [22], respectively.Whereas, the CDC of the weathered oil was obtained from [41].For each layer, the effective CDC was computed using the refractive mixture formula [22].The upper half-space is air with real dielectric constant, ε 0 = 1, while the lower half-space is seawater with CDC, ε w derived from Klein and Swift formulas [22].
In the monostatic setup, the cross-polarization channel within the first-order multilayered SPM is zero; therefore, the results presented in Section III-A include the copolarization channel in VV and HH polarizations.Using only two copolarization channels limits the polarization diversity of the NRCS.Such diversity is crucial for deriving polarimetric parameters used to discriminate scenes with low backscatter contrast or strong volume scattering mechanisms (e.g., see [48] and [49]).However, this limitation is not an issue in this study, which focuses on a bare oil-covered ice scene predominately influenced by surface scattering mechanisms.The cross-polarization backscatter in this case is very low and could be at the NESZ level, making it not usable for inversion.As a result, in an operational setting under our study, having copolarization backscatter measurements across multiple incidence angles is more important for inversion than having measurements from all polarizations at a single incidence angle.
To help interpret the effect of surface roughness on the radar scattering, we calculated the copolarization ratio, ϒ VV/HH , which represents the ratio of NRCS in VV to HH polarization.Within the first-order approximation of SPM, in the case of a single rough interface, this ratio is independent of the roughness parameters and only depends on the ice dielectrics and radar incidence angle [23].In addition, to explore how far the microwave radiation propagates through the oil-covered NI profile, we computed the microwave penetration depth for each layer given by the formula in [21].
2) Inverse Solver: Once the scattering properties have been modeled, it serves as the basis for inverse solvers to estimate the unknown physical parameters from the observed scattering data.In this study, the inverse solver relies on the PSO algorithm (a global optimizer inspired by the social behavior of bird flocks, updating their position and velocity to achieve an optimum search solution [50]).
To estimate the oil slick thickness over the NI profile, the PSO algorithm seeks to minimize the objective function, which is a measure of the estimation error difference, ϵ σ o , based on the inputs between the experimental (or synthetic) and modeled backscatter data [see Fig. 1(b)].Thus, the objective function, f , can be expressed in two different ways as (3) and ( 4), shown at the bottom of the page.
In (3), the objective function, f 1 , uses NRCS data for the minimization, where σ o( * ) (in dB units) represents either the synthetic NRCS when testing the performance of PSO or the experimental NRCS when validating PSO.In (4), the objective function, f 2 , uses the ratio of NRCS in VV to HH polarization data for the minimization, where ϒ ( * ) VV/HH (in dB units) follows a similar definition as the NRCS in (3).The symbol, τ represents three low incidence angles (θ ) of 20 • , 22 • , and 25 • when using synthetic inputs.Conversely, τ represents only one incidence angle, 24.5 • , when using experimental inputs.The notations, ⃗ χ and ⃗ χ ′ represent vectors of the unknown and known variables of the oil-covered NI (see Table I for the definitions of these variables).
The unknown variables include d1 , ε 1 , ε 2 , and ε 3 , which are the retrievable variables to be optimized by PSO for the oil-covered NI profile.This indicates that PSO must initialize seven unknowns for the minimization (note that each dielectric medium, ε is complex, with both real and imaginary parts).Meanwhile, the known variables include d2 , d3 , h 01 , h 12 , h 3w , l 01 , l 12 , and l 3w , are a priori information based on the delineation of certain ice types obtained from SAR satellite imagery.Incorporating these input datasets substantially reduces the number of unknowns and improves the estimation accuracy.Nevertheless, we acknowledge that they may not be entirely precise because, in reality, no operational SAR can directly extract the oil/ice interface roughness (e.g., h 12 , l 12 , h 3w , and l 3w ) and underlying ice thickness ( d2 and d3 ) which are obscured by the surface oil cover [17].This compelled us to test the effect of incorrect ice roughness and thickness parameters on the estimation accuracy in our inversion process.
For instance, we will adjust the experimentally measured value of d3 from 4.5 cm to 2.5 and 8.5 cm, and check the estimation accuracy of the 5 mm oil thickness (see Table V in Section IV-B for the adjustments on other subsurface parameter).Our adjustments were constrained by the following: 1) the SPM region of validity: h < 0.05λ , l < 0.5λ , and h/l < 0.3 [35] and 2) the thickness ranges of frazil and nilas, which are 0.1-1 cm and 2-10 cm, respectively [10].
In this inversion framework, the PSO algorithm was implemented in MATLAB [51], with parameterization values provided in Table II.Note that the lower and upper bound of these unknown variables were estimated from [21], [41], [52], [53].We set the lower bound of oil thickness, d1 slightly above zero to ensure that oil must be present before the algorithm works.

IV. RESULTS
Following the implementation of the EM modeling strategy described in Section III, we present our results in two Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.subsections that correspond to our research questions posed in Section I. First, we delve into the NRCS results derived from the forward solver based on the multilayered SPM.Three sets of modeled backscatter datasets with varying slick thicknesses on ice will be shown.Second, we investigate the retrieval of oil slick thickness from the inverse solver based on our PSO-based parametric inversion process.We will showcase two inversion results: the near-ideal case based on synthetic data to evaluate the suitability of the algorithm and the real case based on experimental scatterometer measurement from the SERF 2020 oil-on-ice experiment.

A. Modeled Backscatter Results
Fig. 3 depicts the angular response of our forward solver, which predicted the C-band backscatter from an oil-covered NI at various oil slick thicknesses (ranging from 1 to 7 mm).Nonetheless, only for thicknesses of 1, 3, and 5 mm are the NRCS values well above the NESZ level (i.e., −45 dB) for incidence angles <35 • .
To avoid poor SNR, we set two thresholds: a mild threshold of 2 dB and a harsh threshold of 6 dB above the NESZ [see Fig. 3(a) and (b)].Employing the mild threshold, the results showed that the predicted NRCS corresponding to 1, 3, and 5 mm slick thicknesses remained above the "NESZ + 2 dB" (i.e., −43 dB) for incidence angles <30 • .Using the harsh threshold, a significant portion of the predicted data exceeded the "NESZ + 6 dB" (i.e., −39 dB) for incidence angles <25 • .The only exception was a few points in the HH polarization for the 5-mm slick layer, which marginally fell below this threshold.Concerning the measured NRCS data that represents only the 5-mm oil-covered NI, both VV and HH exceeded the mild threshold.However, when the harsh threshold was imposed, the VV polarization met the limit, whereas the HH polarization fell 3 dB below it.In summary, our SNR analysis indicates that NRCS data at lower incidence angle (≤25 • ) will be reliable for the subsequent retrieval analysis.
In Fig. 3(a), we observe a downward linear trend in the modeled NRCS of the oil-free NI (0 mm oil thickness) as the incidence angle increases from 20 • to 60 • .This behavior is expected, especially with the VV polarization showing a higher backscatter level and a gentler slope from −12 to −19 dB, in contrast to the HH polarization, which exhibits a lower backscatter level and a steeper slope from −13 to −29 dB.These oil-free NRCS data, both in trend and magnitude, serve as baseline for comparing the NRCS data predicted from the oil-covered NI.When we simulated a 1 mm oil slick layer over the NI surface, the expected NRCS values for VV and HH polarizations decreased slightly by 9 and 8 dB, respectively.However, with 3 and 5 mm slick thicknesses, the predicted copolarization NRCS values experienced a significant drop, averaging 19 and 26 dB, respectively.To validate the accuracy of our forward solver prediction, we compare between the modeled and measured NRCS for oil slick thickness of 5 mm.The modeled VV completely captures the experimental measurement, while the corresponding HH compares well within 2 dB difference.In general, there was a clearly defined contrast once oil is introduced into the NI scattering simulations.In Fig. 3(c), we show a different insight into the predicted backscatter response, which is the angular response of the ratio of modeled NRCS in VV to HH polarization (denoted as ϒ VV/HH ).The oil-free ϒ VV/HH values increase from ∽1 to 10 dB, and upon introducing oil, these values inconsistently decrease by 2, 1.5, and 2 dB for slick layers measuring 1, 3, and 5 mm, respectively.Comparing this with the simulated NRCS data leads to two notable observations.First, both the oil-free and oil-covered ϒ VV/HH data exhibit an increasing trend as the incidence angle ranges from 20 • to 60 • .Second, as oil slick thickness increases, there is no discernible contrast.Table III demonstrates the simulated penetration depth when microwave radiation propagates through each layer of the oil-covered NI profile.Based on the results, we expected two behaviors of microwave penetration.First, when the oil thickness is less than 2 mm, the microwave energy successfully penetrates through the saline-oil layer and oiled-frazil layer until it reaches the sea ice medium.Second, when the oil thickness is 2 mm or greater, the microwave penetration becomes limited within the saline-oil medium, leading to NRCS values closer to or below the NESZ level.Despite this, we noticed that the backscatter response at a 25 • incidence angle remained above the NESZ level (including the mild and harsh thresholds) until the surface oil thickness approached 5 mm [see Fig. 3(b)].

B. Inversion Results
Table IV presents the results of the PSO estimation of oil slick thickness and the related CDCs in different layers of the oil-covered NI profile.Hereafter, the term "True" refers to simulated values, as the actual oil thickness ( d1 ) was simulated, except for d1 = 5 mm, which we directly measured during our experimental work.Note that the elapsed run time took less than 5 min to complete.
With the synthetic NRCS data at three lower incidence angles (20 • , 22 • , and 25 • ), we noticed a consistent pattern where the estimated oil thickness generally increased with the true oil thickness.For instance, when the true oil thickness was 1 mm, the estimated value closely approximated 0.84 mm.Likewise, for oil thicknesses of 3 and 5 mm, the estimated values proportionately increased to 2.62 and 6.23 mm, respectively.While the corresponding retrievable CDCs exhibited random variations as the oil thickness increased, they were held relatively close to their true values within a reasonable margin of error.
With the synthetic ϒ VV/HH at the same incidence angles, we noticed an unacceptable discrepancy between the true oil thickness and the estimated oil thickness.Moreover, the corresponding retrievable CDCs remained unchanged as the oil thickness increased, and their estimated values were held relatively close to their true values.With the experimentally measured NRCS data at a single incidence angle (24.5 • ), we noticed that the estimated slick thickness closely approximated the true value.Similar to the synthetic NRCS, the accompanying retrievable CDCs ranged reasonably close to their true values.
Fig. 4 depicts the estimation error difference and the convergence of the objective function based on our primary estimated variable of interest (i.e., oil thickness).In Fig. 4(a), we compare the error difference between the true and estimated oil thickness values.For the synthetic NRCS, the error difference increases with thicker oil layers; specifically, the 1 and 3 mm oil thicknesses were underestimated by 3% and 8%, respectively, whereas the 5 mm oil thickness was overestimated by 25%.For the experimental NRCS, the error difference was found to be 8% overestimated based on the 5 mm oil thickness.
Fig. 4(b) demonstrates how the objective function progressively minimizes the difference between the modeled and experimental NRCS.With each iteration, the function evaluation consistently decreased, eventually reaching the optimal global minimum ( f 1 = 0.000432) at iteration 62.
Table V presents the retrievable variables of the oil-covered NI profile using incorrect subsurface known ice thickness and roughness data.Our primary estimated variable of interest is the 5 mm oil thickness.For the oiled-frazil layer, we observed that the retrieved oil thickness closely matched the true value (within 0.7 mm average difference).For the sea ice layer, we found that the estimated oil thickness closely aligned with the actual value, with an average difference of 0.4 mm.Regarding the subsurface roughness, we tested two cases: 1) independently adjusting either the h (marked with the * superscript) or l (marked with the * * superscript) and 2) simultaneously adjusting both h and l values (marked with the † superscript).At the oil/ice interface, when we make independent adjustments, the retrieved oil thickness closely matched the true value (within 0.5 mm average difference).However, with the simultaneous adjustments, the estimated Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE IV ESTIMATED VARIABLES OF THE OIL-COVERED NI PROFILE
oil thicknesses were deemed unacceptable.On the contrary, at the ice/water interface, all the estimated oil thickness results closely aligned with the actual value, differing by an average of 0.1 mm, irrespective of the independent or simultaneous adjustments.

V. DISCUSSION
Our primary objective in this manuscript has been to describe how C-band radar data can be used to retrieve oil thickness on top of NI terrains.Two questions are raised to address this objective: what is the expected C-band radar backscatter when a simulated oil slick thickness on the NI surface increases?And how feasible is it to retrieve the slick thickness on NI using experimentally measured C-band radar backscatter?It is clear that the answers to these questions touch upon the issues of modeling backscatter sensitivity, which is a forward scattering problem, and oil spill retrieval, which is an inverse scattering problem.

A. Modeled Backscatter Sensitivity
To address the first question, it is important to distinguish between the backscatter characteristics of oil-free NI and oil-covered NI.Our study presents a novel approach of using a first-order multilayered-SPM forward solver to predict the scattering behavior of stratified bare, undeformed NI; nonetheless, remarkable progress has been achieved with other forward solvers ( [12], [31], and references therein).For instance, Nghiem et al. [12] used a multilayer scattering model that considered both surface and volume contributions to predict the C-band scattering response of the 8 cm bare NI profile, following similar stratification shown in Fig. 2(b).Their modeled results predicted copolarization NRCS, which decreases over the range of incidence angles with VV having a gentler slope and a higher value than HH.In comparison, our multilayer scattering predictions, which only account for surface contributions, agreed with the decreasing trend of their NRCS, except that our NRCS magnitude was 10 dB higher.This discrepancy is due to the differences in experimental setups, ambient environmental conditions and scattering models.
In fact, it is reasonable to assert that modeling the C-band backscatter of bare oil-free NI is relatively well-established, and any remaining unclear scattering mechanisms can be readily identified and are likely to be addressed in the near future.However, when the NI surface is covered with frost flowers, snow, or a combination of both, considerably more work remains to be done, particularly modeling of frost flowercovered NI, which is a common phenomenon during freeze-up seasons.
For oil-covered NI, we performed a backscatter sensitivity study to assess the impact of varying oil layer thickness on bare NI through backscatter simulations.The results in Fig. 3(a) and (b) revealed that the predicted signatures exhibited a highly sensitive response with decreasing NRCS values as the oil layer thickness increased.We attributed this to signal attenuation within the saline-oil layer due to a very large imaginary part of the saline-oil dielectrics.Past experiments by Firoozy et al. [40] and Asihene et al. [17] supported this finding.The study in [40] injected a larger volume of oil (20 L) beneath the ice, resulting in a thicker oil layer (value unknown) when the oil migrated onto the ice surface.Meanwhile the study in [17] injected a smaller oil volume (6 L) beneath the ice, leading to a thinner oil layer (5 mm) when the oil migrated onto the ice surface.Both studies reported a decrease in backscatter sensitivity, with the former demonstrating a more pronounced effect.Although our study is a simulation where oil is spilled on top of ice, the observations from these studies confirm the existence of complex physical processes at interplay between oil and sea ice, which contribute to the smoothening of the NI surface and alteration of its dielectric properties.
For validation purposes, our simulation result compares well with the experimentally backscatter (at 24.5 • incidence angle) of 5 mm oil-covered NI.This was expected as the presence of oil significantly reduces the roughness rmsheight, leading to a smooth oil/ice interface.
Furthermore, we hope to determine how much an increase in oil layer thickness suppresses the backscatter response closer to, and even below the NESZ level (−45 dB) of our C-band microwave scatterometer.This level was met at lower incidence angles (≤35 • ) up to 5 mm oil thickness [see Fig. 3(b)].Despite this, as shown in Fig. 3(a), the backscatter responses of the oil-covered ice (especially the 5-mm thick curve), were in close proximity to the NESZ background, leading to a poor SNR [25].This indicates that a potential noise leakage from the receiver noise floor could compromise the usefulness of NRCS for oil spill detection [24].Using the mild threshold (NESZ + 2 dB), preserve the acceptable signals at incidence angles ≤30 • , while the harsh threshold (NESZ + 6 dB) preserves the useful signals at incidence angles ≤25 • .These thresholds highlight that the NRCS signal values at low incidence angles are well-suitable for analyzing low backscatter phenomena (e.g., [44], [54]).Given our SNR analysis, we recommend the mild threshold for ground-based scatterometer, and the harsh threshold for airborne/spaceborne SAR measurements of oil-covered NI.
Comparing our NESZ level to that of radar satellites with high noise backgrounds (e.g., −33 to −38 dB for RCM ScanSAR low noise beam mode [45]), we speculate that our prediction could be limited to a slick layer of 2 mm or even lower.This highlights the need for technological advancements in the development of new SAR satellites with low NESZ levels, as well as the importance of conducting more in situ radar measurements for calibration purposes.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Other factors we considered include the sensitivity of C-band radar over the range of incidence angles.As noted in Fig. 3(a) and (b), we found that the backscatter sensitivity is greater at lower incidence angles compared to higher incidence angles (>35 • ).This difference in sensitivity arises from the dominance of surface scattering at lower incidence angles, whereas specular reflection becomes more dominant at higher incidence angles.
To establish a relationship between the physical properties and backscatter sensitivity of variations of oil thickness on the dielectric profile of NI, we modeled the penetration depth of microwave propagation at C-band frequency (see Table III).
Our findings indicate that in the absence of a saline-oil cover, microwave radiations can penetrate deeper into the oil-free NI, resulting in a higher backscatter response compared to oil-covered NI.This is primarily due to surface scattering contributions resulting from the strong dielectric contrast at the air/ice interface.On the other hand, when a saline oil is present, its surface weakens the dielectric mismatch and suppresses the surface scattering mechanism from the oil/ice interface.This phenomenon has been described by Asihene et al. [17] as a "masking effect" that leads to a constant backscatter response over time.
Based on our simulated microwave penetration depth data across each layer of the oil-covered NI (see Table III), we initially expected the backscatter sensitivity to approach or even fall below the NESZ level as the oil thickness increases above 2 mm.Surprisingly, we found that the backscatter sensitivity remained above the NESZ level until the surface oil thickness exceeded 5 mm [see Fig. 3(a) and (b)].This indicates that even though the estimated penetration depth through the saline-oil layer is 2 mm (which is defined as the depth where the EM power drops by a factor of e = 2.7 times at zero degree incidence angle [22]), some fraction of EM energy still reaches the oil/ice rough surface at the 5 mm depth.As a result, this phenomenon generates a measurable backscatter signal higher than the NESZ level up to a certain incidence angle.
The sensitivity of the ϒ VV/HH to oil slick thickness variations provides a different perspective on our investigation because of its ability to discriminate between uncontaminated and oil-contaminated NI [43].Prior to the oil spill, ϒ VV/HH exhibited a steep increase over the range of incidence angles.This behavior has been linked to highly saline brine cover under freeze-up conditions, which reflects more microwave energy at larger incidence angles, decreases horizontally polarized wave transmission, and increases the ϒ VV/HH values [55].When oil is introduced, ϒ VV/HH decreases with a clear contrast from the oil-free ice, which we attributed to the reduction of the surface roughness effect.This finding agrees with the observation by Brekke et al. [54] which found that oil reduces the SPM-simulated ϒ VV/HH response.
In comparison to the NRCS, ϒ VV/HH displays inconsistent sensitivity to variations in oil thickness [refer to Fig. 3(a) and (b)].While sensitivity becomes apparent upon oil introduction, it follows neither a decreasing nor increasing trend with further increments in oil thickness.This inconsistency could be attributed to the VV/HH ratio, which is less sensitive to the roughness effect compared to either the VV or HH NRCS.This is the first time, to our knowledge, that such a result has been observed for oil-covered NI; thus, additional experiments data are required for validation.Moreover, when compared to NRCS, the ϒ VV/HH displayed a negligible contrast at lower incidence angles and very low contrast at higher incidence angles as the oil layer thickness increased [see Fig. 3(c)].This indicate that the ϒ VV/HH data was insensitive to a changing oil slick thickness.As a result, we attest that inversion studies involving the retrieval of oil slick thickness on sea ice are more likely to be successful with remotely sensed NRCS observations.

B. Oil Thickness Retrieval
After interpreting the scattering sensitivity predicted by our multilayered-SPM forward solver, we focus our attention on the second question involving the inversion of oil layer thickness on bare, undeformed NI terrains.To tackle this, we used two different backscatter datasets during the PSO-based parametric inversion process to evaluate the data misfit objective function.
Initially, we employed synthetic inputs to test the applicability of our optimizer.The results in Table IV demonstrate that when using NRCS, the estimated oil thickness was closer to the true value compared to when using ϒ VV/HH .This indicates that the PSO hyperparameters were optimally finetuned, and NRCS datasets are suitable for accurately retrieving oil thickness, as their predicted signatures had shown to be sensitive to the increasing oil layer thickness on ice [see Fig. 4(a)].Our findings agree with several previous inversion studies that utilized NRCS datasets for their inversion studies (e.g., [12], [31], [56] and references therein).On the contrary, the inversion using ϒ VV/HH was unsuccessful because the simulation results were inconsistent and lacked contrast with oil thickness variations.We discussed in the preceding Section IV-A that this inconsistency and lack of contrast was caused by the derivation of ϒ VV/HH , which is less sensitive to roughness, making it solely dependent on CDCs of the inhomogeneous medium and radar incidence angle [23].This highlights the significance of roughness parameters for accurate estimation of oil thickness on sea ice; thus, we recommend using instruments such as LiDAR for precise measurement of surface roughness [57], as well as roughness at the underlying oil/ice and ice/water interfaces (although the technology is currently unavailable).
In addition, we used experimentally measured NRCS data to validate the performance of PSO, and the estimated result closely matched the true oil thickness (see Table IV).Compared to synthetic data that employed three incidence angles, the estimation error difference result was impressive with the measured NRCS data at a single incidence angle (see Fig. 4).We acknowledged that using a single dataset in terms of incidence angle and oil thickness is insufficient because our studies were limited to experimental designs constraints.Our future work will aim to improve the realism of the inverse solver by expanding the inversion method to encompass a wide range of incidence angles and oil layer thicknesses.
Our PSO-based parametric inversion technique allows for quantitative extraction of multiple physical parameters from Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.remote sensing observations.Although oil slick thickness is our primary parameter of interest, we were able to retrieve the CDCs of saline oil, oiled frazil, and sea ice (see Table IV).These retrievable CDCs were randomly close to their true values within a reasonable margin of error (except for the 5 mm oil thickness estimation, which we attributed to the noise leakage from the NESZ level).We say that the error margin is "reasonable" because the environmental conditions that contribute to sea ice formation in regions including the Arctic, undergo dynamic processes [58], resulting in constant changes in the thermophysical properties that characterize these inverted CDCs.As a result, we anticipated these random variations, but we recommend that evaluations of the objective function be monitored.For instance, as shown in Fig. 4(b), the convergence of the objective function demonstrates a guaranteed optimal value of our estimated parameter of interest.Moreover, we recommend that future work consider direct dielectric measurements of saline-oil and oil-contaminated ice in a wide range of temperatures and salinities.These measurements could be used for validating the inversion algorithms and also for tuning/validating the dielectric mixture model.
In time-sensitive emergency situations such as oil spills, rapid identification of the thickest oil spill zone is crucial for mounting an effective response.This necessitates the use of fast inverse scattering models capable of extracting the oil thickness from remote sensing observations.Our inversion method takes an iterative approach that uses the forward solver to assess the data misfit objective function, while relying on the PSO algorithm to minimize error differences in the estimated unknown variables.As shown in Table IV, the computational cost of our method is extremely low, taking less than 5 min to complete (on a standard PC with 8 GB RAM and 1.6 GHz 4core CPUs).This run time demonstrates the efficiency of our forward solver and its potential to handle larger radar datasets (such as SAR imagery).
Our inversion studies have been successful thus far; however, it is important to recognize the limitations of this method.First, this inversion process does not support a zerooil thickness.We assumed that oil must be present, which implies that the occurrence of an oil spill situation should be detected before its application.Second, this inversion method is only generalized for oil covering bare NI and does not account for the presence of frost flower and snow cover, which are commonly associated with NI growth [20].This limitation arises from insufficient snowfall and unfavorable atmospheric conditions for frost flower formation during the SERF oilon-ice experiment [17].Our future work will consider this limitation and introduce surface cover such as frost flower, snow, or combination of both into the multilayer dielectric profile NI.Third, this inversion method relies on prior knowledge associated with the inhomogeneous multilayered oil-covered NI (e.g., roughness, ice thickness, and water dielectrics) to improve estimation accuracy.In practice, the existing radar satellites such as SAR, lack the capability to directly retrieve (at least with the limited number of measurements) some of these input parameters, including the roughness at the oil/ice and ice/water interfaces, as well as the underlying ice thicknesses of the oil-covered NI (see Fig. 2).This implies that there is a likelihood of using inaccurate a priori information, which could lead to estimation errors.
For this reason, we tested the estimation accuracy of our inversion with erroneous a priori subsurface roughness and ice thickness (see Table V); note that no testing was done with incorrect surface roughness because oil spills on NI scenes usually smoothens the ice surface [17].Our findings emphasize the importance of using accurate a priori knowledge for precise oil thickness retrievals in oil-covered NI profiles.Understanding these relationships, we establish that PSO-inversion method achieves optimal results when NI thickness is between 2 and 10 cm, with the topmost layer (frazil) thickness falling between 0.1 and 1 cm.Furthermore, we found that the roughness at the oil/ice interface can significantly affect the inversion process compared to the ice/water interface.This is because the presence of highly saline-oil cover attenuates high microwave signals, limiting the microwave propagation through the sea ice medium (see Table III).

VI. CONCLUSION
In this manuscript, we presented an EM modeling strategy for retrieving the thickness of an experimental oil spill on sea ice.Our approach involves using a first-order approximation of a multilayered-SPM forward solver to simulate the C-band radar backscatter from an oil-covered NI and employing a PSO-based parametric inverse solver to estimate the oil slick thickness.
Prior to our simulation and inversion analyses, we checked the influence of noise contamination using mild (2 dB) and harsh (6 dB) thresholds of the measured and simulated NRCS data above the NESZ level (see Fig. 3).Our results showed that NRCS data for the 1, 3, and 5 mm slick layer curves at lower incidence angles ≤25 • were reliable for analysis.This finding highlights the need for technological advancements in the development of new SAR satellites with low NESZ levels.
Our simulation results demonstrate several important findings (see Fig. 3).First, the predicted backscatter is highly sensitive to the increase in oil layer thickness, resulting in a continuous reduction in backscatter strength.This reduction is attributed to signal attenuation within the saline-oil layer.Second, the copolarization ratio is initially sensitive as soon as oil covers the NI surface, decreasing its magnitudes.However, unlike the backscatter sensitivity, the copolarization ratio lacks consistency due to the reduction of the roughness effect.Third, the penetration depth through each layer of the oil-covered NI is significantly hindered by the oil layer, which contains highly saline brine wicking from the topmost layer of NI.
Regarding our inversion results, the estimated oil thickness closely matches the true value (see Table IV).For instance, the synthetic NRCS underestimates 1 and 3 mm oil thicknesses by 3% and 8%, and overestimates 5 mm thickness by 25%, while the experimental NRCS overestimates 5 mm oil thickness by 8%.The corresponding retrievable dielectrics exhibit random variations within a reasonable margin of error, attributed to the dynamic nature of NI thermophysical properties.With the use of synthetic copolarization ratio, the estimated oil thickness proves unsuccessful, and we recommend against its use.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
These findings hold significant implications for oil spill response efforts.Accurately estimating the thickness of an oil slick allows for the identification of areas with the highest concentration of oil, which should be prioritized for immediate action using remediation techniques such as in situ burning and chemical dispersants.By swiftly providing this information, the disaster response team can offer actionable methods and recommendations for containing and cleaning up the spill.In the context of oil spills, time is of the essence, and the rapid identification of the thickest oil spill zones is crucial for making well-informed decisions that minimize damage to natural resources, communities, and economies dependent on them.Additionally, combining these thickness estimates with the extent of the spill area enables more precise calculations of the total volume of oil spilled, facilitating accurate assessments of environmental impacts.
While our study establishes the foundation for using C-band radar observation to estimate oil slick thickness on NI, there is still room for improvement.First, our inversion method relied on a single experimental data set for validation.Second, the method assumes the presence of oil, making it applicable primarily as a planning tool for oil spill situations.Third, the a priori parameters used in our inversion process cannot be directly extracted from current radar satellites.However, this is not a problem provided our assumptions for subsurface NI thickness are within 10 cm, and for underlying roughness remain within the SPM region of validity (see Table V).Finally, surface features such as frost flower and snow, which commonly occur during sea ice formation, were not considered in our inversion method.In future research, we intend to address these limitations by incorporating a more generalized scenario encompassing various surface features and a wider range of experimental data.

Fig. 1 .
Fig. 1.Overview of electromagnetic modeling strategy used for estimating oil thickness on newly formed sea ice: (a) forward solver and (b) inverse solver.

Fig. 3 .
Fig. 3. Modeled C-band backscatter from an oil-covered NI: (a) angular response of co-polarization NRCS at varying oil thickness, (b) thickness response of co-polarization NRCS at 25 • and 50 • incidence angles, and (c) angular response of co-polarization ratio.For comparison, we included the experimentally measured C-band VV and HH polarizations (VV (meas) and HH (meas) ), collected from a 5-mm-thick layer of oil on top of the established NI [17].

Fig. 3 (
b) shows the oil thickness response of the modeled NRCS at incidence angles of 25 • and 50 • .As expected, the NRCS exhibited a continuous decline with increasing oil thickness, and also the predicted VV and HH signatures were more closely related at 25 • compared to the 50 • incidence angle.

Fig. 4 .
Fig. 4. Inverted oil slick thickness: (a) estimation error difference between the true and estimated oil thickness values, based on the synthetic (Syn) and experimental (Exp) NRCS data.The vertical line is a separation between "Syn" and "Exp" data; and (b) objective function evaluation.

TABLE I INPUTS
DATA FOR THE ELECTROMAGNETIC MODELING STRATEGY

TABLE III SIMULATED
MICROWAVE PENETRATION DEPTH THROUGH EACH LAYER OF THE OIL-COVERED NI PROFILE

TABLE V ESTIMATED
VARIABLES OF THE OIL-COVERED NI PROFILE USING ERRONEOUS SUBSURFACE ICE THICKNESS AND ROUGHNESS