Oil Spill Thickness Determination from L-band Synthetic Aperture Radar

The capability of synthetic aperture radar (SAR) to identify oil spills on open water has been known and demonstrated since the early 1970s (Guinard, 1971). SAR is now used routinely to detect accidental and intentional mineral oil releases from vessels, pipelines, and production and extraction facilities in ocean waters. A good example is the National Oceanic and Atmospheric Agency (NOAA) National Environmental Satellite, Data, and Information Service's (NESDIS's) Marine Pollution Products, which report likely spills based on any available optical or radar images (https://www.ospo.noaa.gov/Products/ocean/marinepollution). Similar watchdog services are offered in other countries (e.g., Carpenter, 2016). NOAA's focus is U.S. coastal waters, particularly the Gulf of Mexico, but NOAA also supports observations of spills in international waters when requested. The NOAA Marine Pollution Surveillance Reports show the remote sensing image, identify possible oil and possible thicker oil, provide information about wind conditions, and include notes on the possible source if known.

Despite the need for more extensive information about slick properties, a satellite remote sensing capability for determining oil type or slick thickness and volume remains elusive (Brekke and Jones, 2019). The fact that the capability has not yet been operationalized is primarily due to oil slicks being radar-dark, so that satellite SAR backscatter from slicks is often contaminated by instrument noise to the point that these slick properties cannot be unambiguously determined (Alpers et al., 2017; Angelliaume et al., 2018; Espeseth et al., 2020). Nevertheless, by now a number of studies indicating that thickness or volumetric oil fraction (oil-to-water ratio) are measurable with SAR have been published (see, e.g., Wismann et al., 1998; Minchew et al., 2012; Garcia-Pineda et al., 2013; Boisot et al. 2018; Jones and Holt, 2018). Most studies compare measurements to empirical scattering models such as the Tilted Bragg Model (Valenzuela, 1978) or other models that simplify the ocean surface by considering only a small part of the wave spectrum, and sometimes only a single component. These models ignore the complexity and dynamics of the ocean surface, so are incomplete and introduce uncertainty in the modeling process.

The ocean spectra damping is not linearly related to the radar backscatter intensity, so a full simulation is needed to determine the impact of oil on the measured radar signal. Many evaluations of SAR backscatter from slicks consider only the single ocean wave spectral component corresponding to the Bragg wavenumber at the radar's center frequency. However, this is a simplified assumption that does not account for the complexity of a realistic ocean surface where tilts can readily change the local Bragg scattering component. Also, in addition to slicks modifying the amplitude of the wave spectra, the wind field over a slick is modified and that changes the wave spectra independent of the oil damping factor. A full simulation of electromagnetic scattering based on Maxwell's equations and using a realistic ocean surface can account for backscatter from the full ocean wave spectrum with or without oil damping.

The goal of our study is to determine how accurately slick thickness can be determined from L-band SAR, including the extent to which weather conditions change the SAR backscatter and affect the parameters used to determine the oil thickness. To address the complexity of modeling oil slicks in the open ocean, we developed a physics-based electromagnetic scattering model using a realistic ocean wave spectrum and an established damping function that models the effect of the oil slick on the ocean wave spectrum. L-band SAR data acquired by the Uninhabited Aerial Vehicle SAR (UAVSAR) sensor (Hensley et al., 2009) was used to validate the model for clean ocean and to compare model output to oil slick measurements for a variety of wind conditions. Using the data and the model, we evaluated the extent to which external calibration is needed for the SAR thickness determination, including how changes in wind condition alter the measurements and whether recalibration is needed as the oil in the slick weathers on the ocean surface.

The objectives to meet this goal are (1) determine the wind-dependence of oil-induced ocean wave damping; (2) evaluate the conditions under which relative and absolute thickness can be determined from SAR; (3) develop an algorithm for relative oil thickness determination from SAR data; (4) determine what independent in situ or other measurements are needed to calibrate the parameters of the algorithms; and (5) evaluate the sensitivity of the algorithm to the expected range of ambient conditions and oil properties. We use the Bonn Oil Appearance Codes (Table 1) as a guide to accuracy of the current operational capability for determining layer thickness. The Code specifies a range of thickness of approximately one order of magnitude for each category based on the appearance of the slicked surface to a human observer. Our project is constrained by the limited set of meteorological conditions covered by UAVSAR oil spill and clean ocean acquisitions, and the accuracy evaluation is constrained by the limited availability and accuracy of in situ measurements of thickness acquired concurrent to the UAVSAR acquisitions.

For the work reported here, we evaluate the capability of low noise L-band (1.26 GHz) synthetic aperture radar (SAR) to characterize oil slicks by thickness, both relative and quantitative. The study has both observational and theoretical components. The observational component involved evaluating the existing oil spill data acquired by UAVSAR sensor to determine three things: the range of measured NRCS values under different wind conditions, the contrast in NRCS between slick and water, and the variation of the NRCS within the slick. The theoretical component of the work involves realistic physics-based modeling of the damped and undamped sea surface accounting for all components of the ocean wave spectrum.

The study includes validation campaigns for the developed algorithms using UAVSAR at the Santa Barbara seep field and during a Norwegian oil release and recovery activity, but those campaigns have not yet happened. They are planned for calendar years 2021 and 2022.

Physical models were developed that could account for radar backscatter from seawater with or without an oil layer. Simulated backscatter experiments are run to quantify the effect of slicks composed of oil with different properties and thicknesses. The model flow is shown in Figure 1. Separate components of the model are described below.

The developed electromagnetic scattering models implement Maxwell's equations for 1-D and 2-D simulated ocean surfaces generated using the Elfouhaily ocean wave spectrum (Elfouhaily et al., 1997). An example of the spectrum and simulated sea surface is shown in Figure 2.

We developed 2-dimensional (1-D scattering surface) and 3-dimensional (2-D scattering surface) polarization-dependent electromagnetic models, and cross-calibrated the two models to verify that the 2-D model could be used to accurately model microwave backscatter from seawater with or without an oil layer, and for different thicknesses and dielectric properties of the layer. There are several available electromagnetic scattering models to compute scattering from 1D and 2D rough surfaces. For 2D electromagnetic scattering from 1D sea surfaces, we selected the 2D Finite Difference Time Domain (FDTD) method (Yee, 1966), a previously validated numerical method. The model is based on updating Maxwell's equations on a discretized grid. This model can solve scattering from any dielectric distribution within its computational domain. This provides the flexibility to study scattering from two coherent rough surfaces with different physical properties that are in contact, and thereby quantitatively evaluate the effect of the oil's dielectric properties for a pure oil or oil emulsion layer overlying clean seawater. Though the numerical solution is computationally expensive, the cost of its 2D solver is acceptable for our analysis.

Computing 3D electromagnetic scattering from large 2D ocean surfaces is very challenging for L-band irradiation (23.8 cm wavelength) because scattering at that scale is sensitive to both large-scale and small-scale features of the ocean surface. The asymptotic analytical solutions, e.g., Small Perturbation Method (SPM), Small Slope Approximation (SSA), etc., have validity domains that are dependent on the ocean roughness and radar incident angles. The Kirchhoff Approximation (KA) solver represents sea surfaces as facets. To cover the small-scale features present under conditions of surface damping, the number of small facets can easily become too large to be practically computed. Numerical solutions such as 3D Method of Moments (MoM) or FDTD require extremely large computational resources and become impractical to use for analysis. Therefore, in this work we investigated the 3D Stablized Extended Boundary Condition Method (3D SEBCM) (Duan and Moghaddam, 2012), which is a semi-numerical approach that can cover large roughness range and provide full wave solutions in all radar polarizations. Despite having higher computational efficiency, the 3D SEBCM solver-based Monte Carlo simulation requires much longer run time than the 2D FDTD simulation. Therefore, we used 2D FDTD as the baseline model for sensitivity analysis and developing oil thickness retrieval approach, and performed comparisons between 3D SEBCM and 2D FDTD of a few representative cases to understand the difference between the damping effect in 2D and 3D scattering simulations, and to calculate and apply necessary adjustments when using the 2D FDTD to model the ocean surface.

L-band SAR data from the UAVSAR sensor were used to validate the 1-D surface model against measurements with the Aquarius scatterometer (Yueh et al., 2014) for clean ocean and a variety of wind conditions (Figure 3).

We selected the Jenkins damping function (Jenkins & Jacobs, 1997; Pinel et al., 2010) to model modification to the ocean wave spectra due to the presence of an oil layer because it includes dependence on layer thickness. The Jenkins damping function, also known as the Model of Local Balance, incorporates the wind growth rate effect on the ocean spectra, the effect of oil on wind velocity, interfacial properties at the air-oil interface and the oil-water interface, and oil and seawater bulk properties. Explanations of the damping effect of oil on the sea surface and derivations of the damping functions can be found in (Pinel et al., 2010; Jenkins and Jacobs, 1997; Ermakov et al., 1992; Lombardini et al., 1989). More details on the adjustable parameters and their values for oil slicks are given in Section 4.

The damping ratio, or VV-intensity contrast between clean and slicked water, was selected as the parameter to use for developing an algorithm for determining relative thickness based on previous work that showed the parameter to be very sensitive to oil on seawater (Minchew et al., 2012; Espeseth et al., 2017). Based on the UAVSAR data, a semi-automated method for classifying oil by its relative thickness was developed. The algorithm, which is computationally simple and based on the statistics of the measured values, involves the following steps:

1. Rough segmentation of the scene into oil and clean sea pixels to determine the clean sea NRCS

2. Calculation of the damping ratio (equation (2))

3. Selection of threshold to separate oil from clean sea

4. Calibration against in situ or other data to differentiate sheen from thicker oil classes

For step 1, the probability distribution function (PDF) of the VV NRCS values for the entire scene, including both oil and clean sea pixels, is calculated and the peak associated with clean sea identified as the one with higher values and, for most cases, more pixels. Values around the center of the peak are averaged to determine ⁠. For sensors imaging across a wide range of incidence angles, this calculation is done for bins of ~2° in incidence angle to account the strong incidence angle dependence of the backscatter. This value is then used in Equation 2 to calculate the damping ratio for all pixels in the scene. The PDF of the damping ratio for all the pixels, containing both oil and clean water, is used to identify oil by setting a threshold between the main peak, centered on a value of 1, and the tail of the positive distribution, corresponding to the slick. Alternatively, if multi-polarization data is available, the ratio (co-polarized NRCS ratio) can be used to separate the clean water from oil, again using the peaks in the PDF of this parameter, then the damping ratio used to identify thicker oil as those pixels with higher damping ratio values. The results can be calibrated based on field observations of a slick to identify sheen and other Bonn classes if present, or alternatively without calibration data using a general rule, such as the commonly applied adage that 90% of oil is within 10% of the slick area, to select the thin vs. thick oil threshold. In the event that the measured damping ratios cover a large range (long tail of the PDF), further relative thickness segmentation can be done, as shown in Figure 4 for a large seep in the Mississippi Canyon Block 20 of the Gulf of Mexico imaged with UAVSAR in November 2016. In this case, the identified thicker oil collected along fronts, as one would expect. The method developed is computationally simple, can be applied to single-polarization SAR data, is equally applicable to airborne and spaceborne L-band SAR instruments, and in principle is applicable to frequencies other than L-band. The only step not automated in our implementation is the selection of the threshold in step 3 above, and that can in principle be automated using peak detection and separation algorithms.

A sensitivity analysis was performed using the electromagnetic scattering model to determine which physical properties of the oil and oil-water and oil-air interface most influenced the radar backscatter in the sense that changes in those properties led to the largest changes in backscatter intensity. Using the model, we evaluated whether it is practically feasible to determine the slick thickness in the absolute sense from SAR alone or with a limited set of calibration data.

The Jenkins ocean wave damping model was chosen as the baseline for calculations of the slicked ocean surface because it is physically realistic and models layers of different thickness, including its damping very thin layers equivalent to sheen. Table 2 shows nominal values of the bulk oil and water properties and the interfacial properties at the oil-water and oil-air interfaces. The 2D electromagnetic scattering model was run in a Monte Carlo simulation with different instances of the ocean wave spectra damped according to the Jenkins formula. This was repeated for a physically realistic range of each property listed in Table 2. This allowed us to determine that the surface elasticity (at the air-oil interface), the kinematic viscosity of oil, and the surface viscosity (at the air-oil interface) are the parameters that most influenced the backscatter. Given this information, more simulations were done for the identified parameters to determine the change in backscatter as a function of parameter value, thickness, and incidence angle. The NRCS varies greatly within the range of reasonable values of these parameters, and also varies with the layer thickness. Variations of the NRCS are not linearly related to variations of the oil parameters, and there is no single-valued solution for layer thickness given the range of variation that can occur as the oil in a slick weathers.

We conclude from this study that the relative thickness of oil is represented by higher damping ratio values, as measured with L-band SAR, for typical slick thickness from sheen to thick emulsions. Although not yet validated in field campaigns, the results for relative thickness determined with the semi-automatic damping-ratio-based algorithm presented are physically reasonable. Validation campaigns using the UAVSAR instrument and collections of in situ information are planned for 2021–2022.

Categories of relative thickness can be empirically identified, and calibration could yield Bonn-like classes. A physics-based electromagnetic scattering model indicates that, of the material properties, the damping ratio depends most critically on the surface elasticity at the air-oil interface, the surface viscosity at the air-oil interface, and the oil's kinematic viscosity. The dependency is sufficiently large and complex that given the likely range of values attained through weathering, no one-to-one correspondence between damping ratio and thickness can be made. It is unlikely that the material properties for a slick could be sufficiently well known in an open ocean setting for a rigorous determination of thickness to be made absent a lot of information about the oil and interfacial properties of the slick. However, as Figure 4 shows, relative thickness derived from a semi-automatic process, which could be fully automated, is consistent with what one would expect: sheen near the edges and further from the source, and concentration of oil, possibly emulsified, along fronts and other convergence features.

This information benefits three different activities related to oil spill response, specifically determination of relative thickness to direct skimmers and, assuming the availability of some calibration information, determination of absolute thickness to determine the amount of dispersant to apply, and determination of total release volume given that SAR is already used to determine slick extent. The results can be used for developing a logistical strategy for using low noise SAR to direct and inform spill response.

This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration (NASA). Funding was provided by the U.S. Dept. of Interior Bureau of Safety and Environmental Enforcement (Interagency Agreement E17PG00036) and the NASA Applied Sciences Disasters Program (Task 80NM0018F0830). UAVSAR data are courtesy of the Jet Propulsion Laboratory and can be downloaded from uavsar.jpl.nasa.gov or from the Alaska Satellite Facility (www.asf.alaska.edu).