Internal Waves at the UK Continental Shelf: Automatic Mapping Using the ENVISAT ASAR Sensor

: Oceanic internal waves occur within stratiﬁed water along the boundary between water layers of di ﬀ erent density and are generated when strong tidal currents ﬂow over seabed topography. Their amplitude can exceed 50 m and they transport energy over long distances and cause vertical mixing when the waves break. This study presents the ﬁrst fully automated methodology for the mapping of internal waves using satellite synthetic aperture radar (SAR) data and applies this to explore their spatial and temporal distribution within UK shelf seas. The new algorithm includes enhanced edge detection and spatial processing to target the appearance of these features on satellite images. We acquired and processed over 7000 ENVISAT ASAR scenes covering the UK continental shelf between 2006 and 2012, to automatically generate detailed maps of internal waves. Monthly and annual internal wave climatology maps of the continental shelf were produced showing spatial and temporal variability, which can be used to predict where internal waves have the most impact on the seabed environment and ecology in UK shelf seas. These observations revealed correlations between the temporal patterns of internal waves and the seasons when the continental shelf waters were more stratiﬁed. The maps were validated using well-known seabed topographic features. Concentrations of internal waves were automatically identiﬁed at Wyville-Thomson Ridge in June 2008, at the continental shelf break to the east of Rosemary Bank in January 2010 and in the Faroe-Shetland Channel in June 2011. This new automated methodology has been shown to be robust for mapping internal waves using a large SAR dataset and is recommended for studies in other regions worldwide and for SAR data acquired by other sensors. A.A.K. and P.I.M.; writing—original draft preparation,


Introduction
Oceanic internal waves (IWs) usually develop in stratified water, in which pycnoclines form boundaries that separate water layers of different density. The density difference may be driven by temperature (thermoclines) or salinity (haloclines), with thermoclines occurring more commonly in the open ocean. IWs are common in continental shelf regions, and where brackish water overlies saltwater at the mouths of large rivers. IWs can be generated by strong tidal currents that flow over a seabed of varying topography, such as shelf breaks and shallow sills. IWs are responsible for transferring energy between the large-scale tides and small-scale mixing. On the continental shelf, when their energy reaches the seabed, the shear stress causes sediment disturbance and resuspension. Areas prone to such disturbance form an ideal habitat for certain benthic animals, and hence should be considered for protection when planning offshore developments such as renewable energy devices and oil and gas platforms.
A deeper understanding of the mechanisms of IW formation and interaction with the coastal environment requires the collection and analysis of large volumes of experimental measurements. Traditionally, measurements of IW fields have been carried out using instruments deployed in the ocean, such as moored buoys equipped with temperature, salinity and ocean current sensors or by using acoustic instruments such as sonars on ships. However, in situ measurements are time-consuming, expensive and unsuitable for large-scale observations. Remote sensing (RS) methods provide an alternative to in situ observations of IWs. RS methods are inexpensive and can provide excellent spatial coverage and consistency of observations. They provide insight into the space-time distribution of IWs that are dynamic and diverse.
Synthetic Aperture Radar (SAR) sensors have demonstrated excellent capability for observation of IWs from space. SAR sensors operate by detecting the return of a pulse of radiation in the microwave band and are insensitive to cloud cover and solar illumination. Hence, they can provide all-weather, day and night monitoring of IWs on the ocean surface. Their operation principles are based on high sensitivity to small-scale roughness of the ocean surface induced by wind. The manifestation of IWs in SAR images can be seen as bright and dark bands. Two mechanisms are involved in the generation of intensity bands: modulation of the surface roughness by surface current associated with IWs; and damping of the surface capillary waves by oils and films that tend to concentrate in the regions of surface water convergence, where they form slicks.
Rapid growth in the volume of data produced by satellite SAR sensors over the last decade has stimulated the development of new methods for automatic analysis of SAR images. Several approaches have been developed in the literature for the detection of IWs in SAR data. Rodenas and Garello [1] used a 2-D spectral analysis based on short-time Fourier transform (STFT) to study IW packets in ERS-1 SAR images and to estimate packet wavelength. A similar technique was applied by Changbao et al. [2] to evaluate the dominant wavelength and propagation direction. Rodenas and Garello [3] suggested that the wavelet transform (WT) can perform better than STFT in locating irregularities of the IW signal. They used an analytical IW model to develop a set of wavelet basis functions that were adapted to the IW signatures and applied these functions for analysis of IW patterns in ERS-1 SAR images of the Strait of Gibraltar. The approach based on WT was further extended and a fully automatic tool was developed for detection and orientation estimation of IWs in the ERS-1 SAR images in [3]. The tool employed decomposition of SAR images into different scales using 2-D dyadic WT to suppress speckle noise and improve the accuracy of edge detection. The output of the multi-resolution edge detector was further processed to discriminate look-alike features such as ship wakes, oil slicks and currents. Rodenas and Garello [3] demonstrated the efficiency of this method by processing individual ERS-1 and RADARSAT SAR scenes.
A multiple-stage methodology for IW detection using stationary wavelet transform (SWT) was proposed by Simonin et al. [4], including pre-processing, edge extraction and edge analysis stages. To discriminate look-alikes from the IW features, Simonin et al. [4] implemented an analysis of geometric features of IW edges. However, the algorithm was tested on a very limited set of images and a more detailed study was needed to select IW detection parameters and optimise its performance in various conditions. Bao et al. [5] suggested a convolutional neural network approach for automatic detection of IWs in SAR images produced by the ENVISAT Advanced Synthetic Aperture Radar (ASAR) sensor. The neural network was trained using a database of 950 scenes of IWs in the South China Sea region. The system demonstrated overall good performance but incorrectly identified river estuaries and coastline features [5].
Hence, several advanced approaches have been developed for automatic detection of IWs in satellite images, but none of these have yet been applied to routine processing of large datasets and automatic generation of IW maps. In contrast, a number of SAR studies of the geographical distribution of IWs were based on manual/visual analysis of satellite images. Dokken et al. [6] studied IW activity along the coast of Norway using more than 2600 ERS and RADARSAT-1 images, though the IWs were not combined as a map. Zhao et al. [7]  South China Sea manually traced on ERS-1/2 and Envisat SAR data. Lorenzetti and Dias [8] identified over 460 IW packets in the southeastern continental shelf of Brazil using visual analysis of 264 Envisat SAR images. Kozlov et al. [9] used SAR data to map spatial statistics of short-period IWs in the Kara Sea, based on 89 ENVISAT ASAR images acquired between July and October 2007.
Many studies use a multi-sensor approach for investigation of IWs from space [10][11][12][13][14], based on the joint use of SAR and optical sensor data. Lavrova and Kostianoy [10] processed and analysed SAR and multispectral sensor images to study the temporal and spatial variability of IWs in the Caspian Sea. Several years of Envisat ASAR, Sentinel-1A,1B SAR, Sentinel-2 MSI and Landsat-8 OLI scenes were processed to build a composite map of IW packets to investigate their behaviour with changes in sea surface temperature. Zhou et al. [11] built annual, seasonal and monthly distribution maps of IWs in the Andaman Sea by combining MODIS Terra/Aqua optical and Sentinel-1 SAR images. The SAR sensor was mainly used to fill the gaps in optical data during complex weather conditions. Ning et al. [12] applied data acquired by Sentinel-1 SAR and Gaofen-3 SAR sensors in different periods of time to extend the interval of IW observations in the Malacca Strait. The SAR images were analysed manually to estimate the amplitude and group velocity of IWs. These parameters were linked to the seabed topography in the Malacca Strait. Mitnik et al. [13] applied SAR images acquired by Envisat ASAR, Sentinel-1 and ALOS PALSAR satellite sensors at various polarisations to study IWs and other mesoscale dynamic features generated by the Oyashio current in the northwestern Pacific Ocean. Satellite visible and infrared images were used to complement SAR data. Observations of IWs in the Flores Sea, Indonesia were studied using the Himawari-8 geostationary satellite optical sensor, in combination with MODIS optical sensor images and Sentinel-1 SAR images [14]. Our literature review did not discover any previous studies of the distribution of IWs on the UKCS using satellite data, but studies have been carried out in other parts of the world, where internal waves are frequently observed [6,9,11,12,14].
The review of published literature shows that the study of the geographical distribution and other characteristics of IWs can benefit from the joint processing of satellite data acquired by multiple sensors. Satellite images of IWs were mainly processed manually to minimise detection errors, but this approach also leads to large processing time that only increases with the application of different sensors. In this paper, we attempt to solve this problem by implementing a fully automatic approach that eliminates human errors and enables the processing of large volumes of satellite images. By applying statistical techniques to estimate the geographical distribution of IWs, we reduced the impact of image noise and detection errors. We approach the problem of IW detection in SAR images by considering it as an edge detection task. The problem of edge detection has been extensively studied in the fields of image processing and computer vision [15]. The Canny edge detection algorithm is one of the most successful, providing good accuracy, efficiency and simplicity. We adopted this algorithm for the detection of IW features in SAR data. Based on the Canny algorithm, we have developed a methodology for automatic processing of ENVISAT advanced synthetic aperture radar (ASAR) sensor data and applied this methodology to investigate internal waves on the UK continental shelf (UKCS) in 2006-2012. For this time period, we produced maps of internal wave distribution and climatology. Through this novel approach, we aim to realise the full value of existing and future satellite SAR datasets for systematic observational studies of these important physical oceanographic features.

Study Region
The study was carried out for most of the UKCS, including Dogger Bank in the North Sea, the extension of the UKCS region west of Rockall and parts of the Faroe-Shetland Channel. The continental slope in these areas plays an important role in bringing nutrients to the region and supporting a wide diversity of marine species [16,17]. The outline of the UKCS region used in this study is shown in red in Figure 1, which also shows the bathymetry and seabed topographical features.

ENVISAT ASAR Data
IWs were detected on the UKCS by processing the archive of ASAR sensor data acquired in 2006-2012 by the ENVISAT mission and accessed through the ESA EO Data Gateway Service (https://eogrid.esrin.esa.int/). Data were acquired using both vertical (VV) and horizontal (HH) copolarisation. However, only VV polarisation data were applied to IW detection due to higher sensitivity to variations in sea surface roughness [19].
The ASAR sensor can operate using one of five image acquisition modes, of which the wide swath (WS) mode is most appropriate for IW detection in the UKCS as it achieves the best combination of spatial resolution and coverage. In WS mode the sensor swath is 405 km and the spatial resolution of each pixel is about 150 m × 150 m with a 50% overlap, so the pixel spacing is about 75 m × 75 m. These characteristics are sufficient to acquire more than 1000 images of the UKCS each year and detect most IW features in the area. We examined the total number of ASAR overpasses that were available each month in the sample year 2009. The number of overpasses ranged from 51 in December to 109 in March (see Table 1). There are several reasons why the number of overpasses varied each month: changes of sensor observation timetable and modes; sensor failures and

ENVISAT ASAR Data
IWs were detected on the UKCS by processing the archive of ASAR sensor data acquired in 2006-2012 by the ENVISAT mission and accessed through the ESA EO Data Gateway Service (https://eogrid.esrin.esa.int/). Data were acquired using both vertical (VV) and horizontal (HH) co-polarisation. However, only VV polarisation data were applied to IW detection due to higher sensitivity to variations in sea surface roughness [19].
The ASAR sensor can operate using one of five image acquisition modes, of which the wide swath (WS) mode is most appropriate for IW detection in the UKCS as it achieves the best combination of spatial resolution and coverage. In WS mode the sensor swath is 405 km and the spatial resolution of each pixel is about 150 m × 150 m with a 50% overlap, so the pixel spacing is about 75 m × 75 m. These characteristics are sufficient to acquire more than 1000 images of the UKCS each year and detect most IW features in the area. We examined the total number of ASAR overpasses that were available each month in the sample year 2009. The number of overpasses ranged from 51 in December to 109 in March (see Table 1). There are several reasons why the number of overpasses varied each month: changes of sensor observation timetable and modes; sensor failures and degradation of data quality for some overpasses. The average number of possible overpasses was determined by satellite orbit and the latitude that affects sensor track spacing. The density of the Earth's surface observations is significantly greater for higher latitudes. For the ENVISAT ASAR sensor, the frequency of revisiting the same location at latitude 60 • with incidence angle >5 • is six times per 35 days or every 5.8 days [20]. The ENVISAT ASAR overpasses acquired in WS mode and overlapping with the UKCS in Figure 1 were downloaded and stored for further processing. The results of the preliminary analysis of selected SAR scenes are shown in Figure 2, where the revisit time averaged over the years 2006-2012 was plotted against the month number. It is seen in Figure 2 that a sufficient number of observations is available in the UKCS: on average, any given location was observed every 6.5 days (5.4 times in 35 days) and the maximum interval between observations was 7.9 days in April. The climatology data also show a slight variation in the number of available images in different months (maximum 16% deviation from the mean).
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 22 degradation of data quality for some overpasses. The average number of possible overpasses was determined by satellite orbit and the latitude that affects sensor track spacing. The density of the Earth's surface observations is significantly greater for higher latitudes. For the ENVISAT ASAR sensor, the frequency of revisiting the same location at latitude 60° with incidence angle >5° is six times per 35 days or every 5.8 days [20]. The ENVISAT ASAR overpasses acquired in WS mode and overlapping with the UKCS in Figure 1 were downloaded and stored for further processing. The results of the preliminary analysis of selected SAR scenes are shown in Figure 2, where the revisit time averaged over the years 2006-2012 was plotted against the month number. It is seen in Figure 2 that a sufficient number of observations is available in the UKCS: on average, any given location was observed every 6.5 days (5.4 times in 35 days) and the maximum interval between observations was 7.9 days in April. The climatology data also show a slight variation in the number of available images in different months (maximum 16% deviation from the mean).     Overall, this analysis demonstrates that the number of available ENVISAT ASAR sensor measurements can be sufficient for the estimation of IW monthly, seasonal and annual frequency maps.

Data Processing Chain
The multistage and fully functional processing chain has been developed for automatic processing of the archive of ENVISAT ASAR images. The processing chain is shown in the block diagram in Figure 4 with the arrows indicating the data flow. The chain is implemented using the Python programming language and uses the Sentinel Application Platform (SNAP), the ESA opensource platform for Earth observation processing and analysis. Overall, this analysis demonstrates that the number of available ENVISAT ASAR sensor measurements can be sufficient for the estimation of IW monthly, seasonal and annual frequency maps.

Data Processing Chain
The multistage and fully functional processing chain has been developed for automatic processing of the archive of ENVISAT ASAR images. The processing chain is shown in the block diagram in Figure 4 with the arrows indicating the data flow. The chain is implemented using the Python programming language and uses the Sentinel Application Platform (SNAP), the ESA open-source platform for Earth observation processing and analysis.
At the first data processing stage, the ASAR archive was scanned to select overpasses that geographically overlap with the UKCS outlined in Figure 1. Each of the selected scenes was then processed as follows. First, the image was cropped to include only a section that overlaps with the UKCS. This helps to reduce the size of the processed image and minimise computational expenses. The cropped image was then masked to exclude land from further processing. The land mask was generated using the SRTM3 Shuttle Radar Topography Mission terrain elevation data, with 3 arc seconds spatial resolution in good agreement with that of the ASAR images.
Each ENVISAT ASAR dataset consists of many millions of pixels, hence tiling is necessary. We selected a tile size of 2000 × 2000 pixels. There may be three or more tiles along the swath, depending on how much of the orbit intersects the UKCS. At the first data processing stage, the ASAR archive was scanned to select overpasses that geographically overlap with the UKCS outlined in Figure 1. Each of the selected scenes was then processed as follows. First, the image was cropped to include only a section that overlaps with the UKCS. This helps to reduce the size of the processed image and minimise computational expenses. The cropped image was then masked to exclude land from further processing. The land mask was generated using the SRTM3 Shuttle Radar Topography Mission terrain elevation data, with 3 arc seconds spatial resolution in good agreement with that of the ASAR images.
Each ENVISAT ASAR dataset consists of many millions of pixels, hence tiling is necessary. We selected a tile size of 2000 × 2000 pixels. There may be three or more tiles along the swath, depending on how much of the orbit intersects the UKCS.
The next stage in the block diagram in Figure 4 was to extract the contours of IWs in the image, then apply tests to discriminate true IW features from non-IW edges. This was applied to the original level 1 ASAR data to avoid any loss of data due to image resampling and mapping. The intermediate results of IW detection were saved in the PNG image format at native ASAR resolution. The results of IW detection were applied at the next processing stage to produce the IW occurrence maps that show the frequency of IW observation for different regions of the UKCS (Section 4). Any data outside the UKCS were masked out and the images were mapped into the geographic Lat/Lon projection, WGS 84 geographic datum. Maps with pixel size 10 × 10 arc minutes, or 19 km in latitude and 10 km in longitude at UK latitudes, were generated and saved in PNG, GeoTIFF and NetCDF data format. The stages in the processing chain are detailed in the following sections.

Implementing Canny Edge Detection
The Canny edge detection algorithm [15] consists of multiple stages: reducing noise using a Gaussian filter, then calculating the image gradient magnitude and direction. The edges are thinned by suppressing non-maximal gradients, then a final hysteresis stage retains weaker gradients only if they are connected to stronger gradient pixels along the edge.
This classic algorithm has been extended with additional stages for the IW edge analysis. At the first stage, the image background intensity level is estimated by applying a running window median filter of round shape and 60 pixels in diameter to the ASAR image, which is shown in Figure 5a. The filter outputs a smoothed image with removed ridge features and details, as illustrated in Figure 5b. The next stage in the block diagram in Figure 4 was to extract the contours of IWs in the image, then apply tests to discriminate true IW features from non-IW edges. This was applied to the original level 1 ASAR data to avoid any loss of data due to image resampling and mapping. The intermediate results of IW detection were saved in the PNG image format at native ASAR resolution. The results of IW detection were applied at the next processing stage to produce the IW occurrence maps that show the frequency of IW observation for different regions of the UKCS (Section 4). Any data outside the UKCS were masked out and the images were mapped into the geographic Lat/Lon projection, WGS 84 geographic datum. Maps with pixel size 10 × 10 arc minutes, or 19 km in latitude and 10 km in longitude at UK latitudes, were generated and saved in PNG, GeoTIFF and NetCDF data format. The stages in the processing chain are detailed in the following sections.

Implementing Canny Edge Detection
The Canny edge detection algorithm [15] consists of multiple stages: reducing noise using a Gaussian filter, then calculating the image gradient magnitude and direction. The edges are thinned by suppressing non-maximal gradients, then a final hysteresis stage retains weaker gradients only if they are connected to stronger gradient pixels along the edge.
This classic algorithm has been extended with additional stages for the IW edge analysis. At the first stage, the image background intensity level is estimated by applying a running window median filter of round shape and 60 pixels in diameter to the ASAR image, which is shown in Figure 5a. The filter outputs a smoothed image with removed ridge features and details, as illustrated in Figure 5b.
The smoothed image serves as the baseline for estimation of the image noise level and for the selection of the upper and lower thresholds of the Canny edge detection algorithm. The noise level is estimated by applying a multiplicative model where the noise level is proportional to the image intensity [21]. The multiplicative nature of noise in ASAR images is known as the speckle phenomenon, where the coherent combination of reflected signals with different phases and amplitudes produce random fluctuations [22]. The variance of the speckle noise σ 2 sp in ASAR image is defined by the equivalent number of looks L and the average intensity I [23] as: It follows from (1) that the intensity of speckle noise increases with image intensity, so the threshold for Canny edge detection was selected accordingly to match these changes. The upper threshold T U was calculated for every pixel location (r, c) as: where C U is a constant that defines the sensitivity of the edge detection algorithm and its false alarm rate. σ sp (r, c) is the variance of speckle noise at location (r, c). A similar approach has been applied for the selection of the lower threshold T L : The C U and C L constants in (2) and (3) were selected by processing test ASAR images and visual assessment of detection results.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 22 The smoothed image serves as the baseline for estimation of the image noise level and for the selection of the upper and lower thresholds of the Canny edge detection algorithm. The noise level is estimated by applying a multiplicative model where the noise level is proportional to the image intensity [21]. The multiplicative nature of noise in ASAR images is known as the speckle phenomenon, where the coherent combination of reflected signals with different phases and amplitudes produce random fluctuations [22]. The variance of the speckle noise 2 in ASAR image is defined by the equivalent number of looks L and the average intensity I [23] as: It follows from (1) that the intensity of speckle noise increases with image intensity, so the threshold for Canny edge detection was selected accordingly to match these changes. The upper

Algorithm Tuning
The algorithm for IW detection (Section 3.2.1) uses several parameters that affect its sensitivity and the false alarm rate. To achieve better performance and improve the accuracy of detection results, this algorithm has been additionally tuned using test ASAR images of IWs.
The test data we used were the scene in Figure 5a (Celtic Sea, 13 July 2011) chosen to have as many IWs as possible, and all scenes from 14 consecutive days during which no genuine IWs were observed visually (1-14 July 2009). Images of the original data superimposed with detected edges were manually inspected for errors, both false negative (where an IW-like edge was visible) and false positive (where no IW-like edge was visible, or it was not part of a wave train). Images were displayed at each stage of the algorithm, allowing us to see the effect of each stage. Before tuning, the default parameters were: lower Canny threshold C L = 0.02; upper Canny threshold C U = 0.2; smoothing sigma σ = 3.
Based on previous experience of front detection, we tried upper Canny thresholds from 0.3 to 0.9 and a lower threshold of 0.1. The upper threshold of 0.5 detected some clear IW-like features not detected by higher thresholds, whereas 0.3 detected more IWs but significantly increased the number of false detections. We decided to err on the side of ensuring the detection of IWs with more false positives and filter out the false positives at a later stage of processing (Section 3.2.3). Further reducing the lower Canny threshold resulted in many more false positives with no apparent gain in IW detection. We did not try changing the smoothing sigma, as the chosen value was considered to be reasonable. Final parameters after tuning were: lower Canny threshold C L = 0.1; upper Canny threshold C U = 0.3; smoothing sigma σ = 3. The resulting detected edges are shown in Figure 5c,d.

Filtering and Grouping of IWs
This proceeds in three distinct stages following Canny edge detection, as illustrated in Figure 6. The overarching idea behind this process is to join edges that we consider to belong to the same IW train and split edges at sharp discontinuities indicative either of a junction between two IWs or of a non-IW feature. Once this is done, we can select just the features most clearly belonging to IW trains.  In the first stage, Canny-detected edges are divided into clusters of connected line segments (Figure 5d), and clusters containing less than ℎ pixels are discarded. We have found that some scenes contain many very short, isolated fragments that can slow down processing considerably while contributing little or nothing to algorithm performance. The orientation of such fragments is difficult to determine, so they are often associated wrongly with other segments in the third stage below. Line segments are formed by 'walking' along edges, creating a break between segments if the line bends by ≥45° or we arrive at a previously encountered pixel. At triple points (junctions of three or more edges), the two most similarly oriented edges are considered part of the same segment unless they deviate by ≥45°. Each segment is then simplified by representing it as a number of straight facets, with each facet made as long as possible without deviating from the segment by more than 1 pixel. Facets within a segment are then combined into geometric elements (straight sections and circular arcs, depending on which best fits the segment data). Straight sections are drawn from end to end of the combined data, while arcs are drawn between the two ends and an intermediate vertex. If the In the first stage, Canny-detected edges are divided into clusters of connected line segments (Figure 5d), and clusters containing less than N short pixels are discarded. We have found that some scenes contain many very short, isolated fragments that can slow down processing considerably while contributing little or nothing to algorithm performance. The orientation of such fragments is difficult to determine, so they are often associated wrongly with other segments in the third stage below. Line segments are formed by 'walking' along edges, creating a break between segments if the line bends by ≥45 • or we arrive at a previously encountered pixel. At triple points (junctions of three or more edges), the two most similarly oriented edges are considered part of the same segment unless they deviate by ≥45 • . Each segment is then simplified by representing it as a number of straight facets, with each facet made as long as possible without deviating from the segment by more than 1 pixel. Facets within a segment are then combined into geometric elements (straight sections and circular arcs, depending on which best fits the segment data). Straight sections are drawn from end to end of the combined data, while arcs are drawn between the two ends and an intermediate vertex. If the angle between two adjacent elements is ≥45 • the segment is split, and if an arc has a radius less than rmin, it is considered as a discontinuity and the segment is split before and after the arc. Results presented use parameter settings N short = 5 pixels and r min = 5 pixels (about 375 m).
In the second stage, segments that are longer than N short pixels and either within the same cluster or in a neighbouring cluster separated by less than d max pixels are considered for grouping, which takes place if the combination of two elements (one from each segment) into a single straight section or arc results in a better overall fit. Straight sections are drawn from either end of one element to either end of the other, while arcs use any three of the four ends. After connection, groups that are shorter than N long pixels are discarded. The algorithm uses parameter settings d max = 20 pixels and N long = 25 pixels.
In the third stage, parallel groups are connected. For each group A, all other groups B with facet vertices within d seg pixels of any part of group A are found using a KD-tree algorithm [24], then for each group B, the N pixels within d seg pixels of group A are considered. Their median distance d med and the number N close of distances within x max ·d med of d med are calculated. If N close is greater than N min and N close /N is greater than p min , groups A and B are connected. The purpose of this is to find pairs of groups that have long sections separated by a consistent distance. The greater the separation of the groups, the more consistent their distance must be. Groups are only considered to be IWs if they are connected to at least one other group. Figure 7 shows the results of parallel groups identification. The following parameters were used: d seg = 20 pixels, x max = 15, N min = 25 and p min = 0.5.

Distribution of IWs
The factors that adversely affect the accuracy of IW contour maps such as Figure 7 are speckle noise in ASAR images, low contrast of IWs features and possibly false alarms at the output of the edge detection algorithm. To reduce these errors and improve the accuracy of IW mapping, we adopted a gridding approach based on combining the information in space and time. The satellite measurements were mapped into an evenly spaced grid of lower resolution and averaged over a time interval.
At the output of the detection algorithm, described in detail in Section 3.2, the contours of IWs were produced in vector format. To combine vector data spatially or in time, the contours were transformed into a raster format as follows. For each radar image, a raster image of exactly the same size was created filled with zero values, then IW contours (see Figure 7) were plotted on this image by using Bresenham's line drawing algorithm, filling the grid with non-zero values in the locations where IWs were detected. These raster images were then composited in space and time to produce monthly, annual and seasonal composite maps of the distribution of IW occurrence. For this study, we adopted a spatial resolution of 10 × 10 arc minutes, or about 9 × 18 km at 60 • N.
These maps show how frequently IWs were found in radar images. For each grid cell in the map, the frequency of IW occurrence was estimated as the ratio of detected IW pixels to the total number of radar image pixels. The frequency values were calculated for multiple radar images and then averaged in space and time as: where k is the ASAR image index and K is the total number of composed images. i and j are row and column numbers of a grid element in the IW distribution map, and I,J are the total numbers of grid elements. (x,y) are the coordinates of a pixel in the ASAR image and X i ,Y j is the subset of all ASAR pixels contained in grid element (i,j). The total number of pixels in the subset is denoted by N. The value of w x,y is set to one if the IW was detected in the ASAR image pixel with coordinates (x,y), otherwise, it is set to zero.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 22 identification. The following parameters were used: = 20 pixels, = 15, = 25 and = 0.5. Figure 7. Parallel groups of internal wave features. Each set of parallel groups has its own colour. The total number of candidates for IW contours has reduced significantly after filtering and grouping stages.

Distribution of IWs
The factors that adversely affect the accuracy of IW contour maps such as Figure 7 are speckle noise in ASAR images, low contrast of IWs features and possibly false alarms at the output of the edge detection algorithm. To reduce these errors and improve the accuracy of IW mapping, we adopted a gridding approach based on combining the information in space and time. The satellite measurements were mapped into an evenly spaced grid of lower resolution and averaged over a time interval.
At the output of the detection algorithm, described in detail in Section 3.2, the contours of IWs were produced in vector format. To combine vector data spatially or in time, the contours were transformed into a raster format as follows. For each radar image, a raster image of exactly the same size was created filled with zero values, then IW contours (see Figure 7) were plotted on this image by using Bresenham's line drawing algorithm, filling the grid with non-zero values in the locations where IWs were detected. These raster images were then composited in space and time to produce monthly, annual and seasonal composite maps of the distribution of IW occurrence. For this study, we adopted a spatial resolution of 10 × 10 arc minutes, or about 9 × 18 km at 60°N . The total number of ASAR images combined in Equation (4) is different for monthly, annual and seasonal IW distribution maps. For a one-month time interval, the average number of images is K ≈ 6 ( Figure 2). The number of ASAR image pixels contained in each grid element of the IW distribution map depends on the map resolution and coordinates. With the 10 × 10 arc minute resolution adopted here, each grid cell covers around 250 × 130 ≈ 33,000 pixels. Hence, on average, ∼200,000 ASAR pixels can be composed to estimate the frequency parameter F i,j using Equation (4) for each grid element of a monthly IW distribution map. This makes the F i,j estimate more resistant to the false alarms generated by the IW detection technique, described in Section 3.2.

Distribution of IWs
The distribution of IWs in the UKCS was evaluated at different time scales and intervals. Monthly maps of IW occurrence were generated for each month in the years 2006-2012. The monthly maps provide the best time resolution, appropriate for observing ongoing changes in IW distribution.
However, the number of sensor measurements averaged over one month is relatively small (see Section 2.2) and this affects the accuracy of the IW occurrence maps. Selected examples of monthly IW occurrence maps are shown in Figure 8. The frequency of IW in the maps is presented using a colour palette in the range from 0 (white) to 0.005 (dark red). Masked areas are shown in grey. As can be expected, a higher frequency of IW occurrence is observed along seabed topographic features, such as shelf breaks, ridges and troughs (see Figure 1). Figure 8a-    Annual composite maps of IW occurrence are presented in Figure 10 for the years 2006-2012. It is seen from these maps that over the whole year, IWs are mostly present in the Hatton-Rockall Basin and Hatton Bank. The most significant years in terms of IW activity were 2006, 2010 and 2011.
Monthly climatology maps of IW occurrence are presented in Figure 11. The maps were generated by combining the same month over the years 2006 to 2012. These show that a higher intensity of IWs is observed in the spring and summer periods. These observations agree with the seasons when the UKCS waters are more stratified [25]. Monthly climatology maps of IW occurrence are presented in Figure 11. The maps were generated by combining the same month over the years 2006 to 2012. These show that a higher intensity of IWs is observed in the spring and summer periods. These observations agree with the seasons when the UKCS waters are more stratified [25].
IWs can travel hundreds of km from their source, so IW occurrence is not expected to be directly related to bathymetry, but it is of interest to see whether IWs are detected more often at certain depths.

Evaluation Results
We visually inspected all monthly IW occurrence maps to identify patterns of IW occurrence (see Supplementary Material to access the maps). IW-like features are detected consistently throughout the year very close to the coast, though many of these near-coast detections are likely to be false positives due to the interaction of surface waves with the coastal bathymetry. Further from the coast, we can be more confident in our results. There follows a list of the most prominent non-coastal IW detection features seen in Figures 9-11 in the UKCS, from west to east. The features are denoted with letters and their locations are shown in red in Figure 9. Selected features are shown in closeup in Figure 13. Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 22 Figure 11. Monthly climatology maps of internal wave occurrence: (a-l) from January to December.
IWs can travel hundreds of km from their source, so IW occurrence is not expected to be directly related to bathymetry, but it is of interest to see whether IWs are detected more often at certain depths. Figure 12 shows a bar chart of the mean IW occurrence as a function of bottom depth, binned at 20 m intervals. This indicates that IW occurrence shows no depth trend between about 900 m and

Evaluation Results
We visually inspected all monthly IW occurrence maps to identify patterns of IW occurrence (see Supplementary Material to access the maps). IW-like features are detected consistently throughout the year very close to the coast, though many of these near-coast detections are likely to be false positives due to the interaction of surface waves with the coastal bathymetry. Further from the coast, we can be more confident in our results. There follows a list of the most prominent noncoastal IW detection features seen in Figures 9-11 in the UKCS, from west to east. The features are denoted with letters and their locations are shown in red in Figure 9. Selected features are shown in closeup in Figure 13.   in all years and most months, most strongly from March to July. At this depth, and with this consistency of appearance, it is possible that the feature is due to seabed effects on surface waves.
The features identified at the UK continental shelf using IW occurrence maps were compared to observations presented in other studies. The Atlas of Internal Solitary-like Waves [26] indicates IW activity in the Faroe-Shetland Channel and at Rockall Bank and the north and west coasts of Scotland. These locations are in good agreement with features K, D and J, respectively in Figure 9. The spatial resolution of IW occurrence maps in the atlas is at least 10 times lower than the resolution of the maps generated in this study and many features observed here were missing in the atlas. More detailed information was provided through in situ measurements and modelling experiments carried out in [27][28][29][30][31]. Three main areas were identified in the literature as IW hotspots: the Malin Shelf, the Wyville-Thomson Ridge and Jones Bank in the Celtic Shelf.
Detailed studies of IWs at the Malin Shelf were carried out by Small et al. in [27] using acoustic measurements and SAR data. It was shown that in summer IWs are a common feature in this region. These observations also match with the IW hotspot F in this region in Figure 13d. This feature can also be clearly seen in the climatology maps in July and August (Figure 11g,h). Vlasenko et al. [28] carried out numerical modelling of IW at the Malin Shelf and provided further insight into the behaviour of IWs in this region. This study describes two wave processes developing in this area: a tidal beam generated at the shelf break and the bottom trapped internal waves generated by the tidal flow over local seabed features.
Another IW hotspot studied in the literature is the Wyville-Thomson Ridge. The numerical simulation results performed by Hall et al. [29] indicated the presence of internal tides in this region. Vlasenko et al. [30] carried out a detailed analysis of water circulation at the ridge that demonstrated favourable conditions for IWs.
Internal waves in the Celtic Sea have been studied in situ and with numerical modelling. Vlasenko et al.
[31] carried out numerical modelling of the interaction of stratified tidal flows with the seabed at Jones Bank and found that tidal currents can generate two types of IWs. Palmer et al.
[32] carried out in situ observations at Jones Bank and revealed IWs of up to 40 metres in amplitude occurring during spring periods. Feature G in Figure 9 shows a high occurrence of IWs around Jones Bank.

Conclusions
In this work, a methodology has been developed for processing satellite SAR sensor data to map internal waves in the UK Continental Shelf. We have made significant advances on previous approaches by implementing automated detection of the IW features. This has the advantage of speeding up the analysis, enabling thousands of SAR images to be processed rather than tens. It also avoids the tedious and subjective manual/visual detection of IW packets, a critical factor in the generation of an accurate and comprehensive distribution map.
We provide the first detailed and systematic distribution maps of IW occurrence across the UKCS by processing a 7-year time series of ENVISAT ASAR sensor data, acquired by the European Space Agency. Overall, the IW occurrence maps demonstrate good agreement with previous findings and with prominent non-coastal seabed topographic features. The comparison of annual and monthly IW maps reveals high interannual and seasonal variability of IWs in the UKCS.
Features revealed by this work could be investigated further by construction of time series of IW occurrence in a region of interest, followed by manual examination of the scenes contributing to each feature to check for artefacts due to meteorological events (e.g., convective storms or atmospheric gravity waves), or consistent features indicative of a permanent structure in the water. As structures such as wind farms or oil platforms are placed in or removed from UKCS waters, time series analysis may also reveal the effect of these structures on IW occurrence. Future work will address adapting the methodology to exploit SAR datasets produced by other sensors. For example, ESA's Sentinel-1A and 1B paired sensors would offer an improvement in revisit time, spatial resolution and continuity of IW observations. The advance we have demonstrated in this automated and robust algorithm should be widely applicable to studies of the distribution and impact of internal waves in other regions around the world.