Ocean Front Detection with Glider and Satellite-Derived SST Data in the Southern California Current System

This study proposes a method to detect ocean fronts from in situ temperature and density glider measurements. This method is applied to data collected along the CalCOFI Line 90, south of the California Current System (CCS), over the 2006–2013 period. It is based on image-processing techniques commonly applied to sea surface temperature (SST) satellite data. Front detection results using glider data are consistent with those obtained in other studies carried out in the CCS. SST images of the Multi-scale Ultra-high Resolution (MUR) dataset were also used to compare the probability of occurrence or front frequency (FF) obtained with the two datasets. Glider and MUR temperatures are highly correlated. Along Line 90, frontal frequency exhibited the same maxima near the transition zone (~130 km offshore) as derived from MUR and glider datasets. However, marked differences were found in the bimonthly FF probability with high (low) front frequency in spring-summer for glider (MUR) data. Methodological differences explaining these contrasting results are investigated. Thermohaline-compensated fronts are more abundant towards the oceanic zone, although most fronts are detected using both temperature and density criteria, indicating a significant contribution of temperature to density in this region.


Introduction
Ocean fronts are elongated structures characterized by abrupt spatial changes in physical, chemical, and biological properties that mark the boundary between different water masses [1]. These dynamic structures are susceptible to frontogenesis by larger-scale deformation flows [2] and baroclinic instabilities occurring over a wide range of spatial and temporal scales [3,4]. Fronts are important because they are associated with particleaccumulating convergence zones that form rich areas of high biological activity [5][6][7].
Remote observations from space have fostered the development of front detection algorithms from sea surface temperature (SST) and ocean color satellite imagery and have significantly contributed to the study of frontal variability. Various techniques have been used for ocean front detection from SST images. Classic methods include algorithms (Sobel, Roberts, or Laplacian) for detecting horizontal gradient intensity that use the first-or second-order derivative as a criterion to determine discontinuities [8][9][10]. Another technique for detecting the boundary between two water masses is based on the analyses of the frequency distribution (histogram) of an SST image [11]. This method has been improved and used for mapping oceanic fronts in different regions [12][13][14][15][16]. More sophisticated methods combine different techniques; for instance, Shimada et al. [17] identified fronts in great detail near the coasts of Japan using an approach based on the Jensen-Shannon divergence and morphological operations applied to SST images to extract curvilinear frontal features.
The density of frontal structures in Eastern Boundary Upwelling Systems (EBUS) such as the California Current System (CCS) presents temporal variations. Front frequency, used here as the probability of classifying a pixel as a frontal feature [15], has been related to coastal upwelling seasonal variations in the CCS, attaining a maximum FF in summer and a minimum in winter [18]. Statistical analyses carried out by Kahru et al. [15] showed that FF in monthly and seasonal SST images was positively correlated with colder SST and increasing near-surface Chlorophyll-a concentrations, suggesting an increased injection of nutrients into the photic layer. However, their interannual variations in some oceanic regions of the CCS were not always correlated.
Remote measurements of the ocean surface and the development of suitable detection techniques have been used to improve our understanding of the spatio-temporal (2D) distribution of ocean fronts at local (coastal), regional, and basin scales. However, a thorough study of these structures also requires in situ observations, supporting an unambiguous interpretation of the detected fronts with remote satellite measurements. High-resolution in situ thermohaline and optical (i.e., chlorophyll fluorescence) measurements can be provided by autonomous underwater vehicles (gliders), which have been demonstrated to be more convenient platforms for long-term regional sampling than traditional research vessels [19]. Observations from the California Underwater Glider Network (CUGN; spraydata.ucsd.edu/projects/CUGN, accesed on 22 May 2019) have been used for over 12 years to explore the persistent oceanographic features and physical-biological interactions in the CCS at resolutions not achieved before. The observational climatology revealed a wind-driven upwelling front located about 200 km off the Southern California coast, with an associated equatorward current [20]. Inshore, currents flow poleward associated with the recirculating Southern California Bight eddy [21]. Additional offshore secondary density fronts are abundant in the CCS throughout the year (especially in spring-summer), and they mark abrupt changes in chlorophyll fluorescence and in primary and secondary production [22][23][24][25]. Long-term series from CUGN are also very convenient to address climatological phenomena such as El Niño, which directly affect the variability, location, and extension of the CCS upwelling front and secondary fronts [26].
Continuous measurements recorded by gliders can produce hundreds of vertical profiles of essential oceanographic variables. This makes the identification of frontal structures somewhat laborious if we intend to consider long trajectories (>500 km) such as the Line 90 of California Cooperative Oceanic Fisheries Investigations (CalCOFI), which has been studied since 2006 [27]. To facilitate the identification of frontal zones, this paper introduces an automated method to detect ocean fronts using repeated measurements spanning over seven years and recorded by autonomous underwater vehicles along a CalCOFI transect located south of the CCS. The front-detection method is based on imageprocessing techniques that have been used for front detection using SST data. The fronts detected using glider data were compared with those identified using moderately highresolution SST images. Frontal structures of temperature and density were also identified with glider data and an estimate of the amount of thermohaline compensated fronts.

Glider Data
The data used in this work are a subset of the historical datasets gathered by the CUGN [27] with underwater gliders (Spray glider, https://spray.ucsd.edu, accesed on 22 May 2019). These autonomous vehicles employ variable buoyancy for diving and wings to glide horizontally and perform a sawtooth trajectory through the water column [27,28]. CUGN Spray gliders are equipped with a SeaBird 41CP pumped Conductivity-Remote Sens. 2021, 13, 5032 3 of 18 Temperature-Depth (CTD) sensor to measure pressure, temperature, and conductivity (salinity). Potential density (σ θ ) is computed from absolute salinity (S A ) and conservative temperature (θ) using the Thermodynamic Equation of Seawater-2010 (TEOS-10) MATLAB (https://www.mathworks.com, accessed on 23 August 2021) package.
Spray gliders can perform multiple dives from the surface to a depth of 500 m over several months, completing a cycle in 3 h while traveling less than 3 km horizontally during each dive [20]. We used a total of 70 backward and forward trajectories along the CalCOFI Line 90 (Figure 1) The data used in this work are a subset of the historical datasets gathered by the CUGN [27] with underwater gliders (Spray glider, https://spray.ucsd.edu, accesed on 22 May 2019). These autonomous vehicles employ variable buoyancy for diving and wings to glide horizontally and perform a sawtooth trajectory through the water column [27,28]. CUGN Spray gliders are equipped with a SeaBird 41CP pumped Conductivity-Temperature-Depth (CTD) sensor to measure pressure, temperature, and conductivity (salinity). Potential density (σθ) is computed from absolute salinity (SA) and conservative temperature (θ) using the Thermodynamic Equation of Seawater-2010 (TEOS-10) MATLAB (https://www.mathworks.com) package.
Spray gliders can perform multiple dives from the surface to a depth of 500 m over several months, completing a cycle in 3 h while traveling less than 3 km horizontally during each dive [20]. We used a total of 70 backward and forward trajectories along the CalCOFI Line 90 (Figure 1), conducted continuously from October 2006 to December 2013. Line 90 extends from ~25 km to ~500 km offshore, that is, between locations P1 (122.59°W, 31.12°N) and P2 (117.95°W, 33.33°N). The resolution of the glider measurements is a function of the sensor sampling frequency and glider speed. The vertical profiles are vertically averaged (grouped) at 10 m intervals from the surface to 500 m depth. However, we only used conservative temperature and potential density data from the upper 10 m for comparison with SST satellite data. Horizontally equidistant along-track measurements were obtained using a third-order polynomial function fitted to irregularly distributed glider locations along each trajectory Equation (1). The longitudes (xi) corresponding to evenly spaced (3 km) locations between P1 and P2 were computed as dx = 3|P1(1)−P2(1)|/3 dist(P1,P2). These coordinates (xi) are replaced in the fitted polynomial function to obtain the adjusted latitudes (yi). The best-fit function obtained was a cubic polynomial of the following form: The resolution of the glider measurements is a function of the sensor sampling frequency and glider speed. The vertical profiles are vertically averaged (grouped) at 10 m intervals from the surface to 500 m depth. However, we only used conservative temperature and potential density data from the upper 10 m for comparison with SST satellite data. Horizontally equidistant along-track measurements were obtained using a third-order polynomial function fitted to irregularly distributed glider locations along each trajectory Equation (1). The longitudes (x i ) corresponding to evenly spaced (3 km) locations between P 1 and P 2 were computed as dx = 3|P 1 (1) − P 2 (1)|/3 dist(P 1 ,P 2 ). These coordinates (x i ) are replaced in the fitted polynomial function to obtain the adjusted latitudes (y i ). The best-fit function obtained was a cubic polynomial of the following form: where x i (i = 1, 2, . . . , N) are the longitudes of the N positions along the trajectory and a 0 , a 1 , a 2, a 3 are the fitted coefficients for a given trajectory. Equation (2) calculates the shortest distance between the glider position and the straight line joining P 1 and P 2 . Once the coordinates are obtained for the f (x i ) = 1, 2, 3... M points, the measured SST value closest to the line (i.e., the shortest Euclidean Distance, D E ) is obtained as follows: Remote Sens. 2021, 13, 5032 4 of 18 which measures the distance to the adjusted points (x a ,y a ) in km. Finally, data values for each trajectory were organized in a matrix containing temperature and density values, distance, and date.

MUR Observations
Global SST images from the Multi-scale Ultra-high Resolution (MUR) data product (http://mur.jpl.nasa.gov, accessed on 30 May 2019) were used for this study to compare fronts identified using in situ glider data and fronts detected from remote satellite observations [29]. These daily, cloud-free, 1-km spatial resolution images combine measurements from various satellites (microwave and infrared) and in situ observations. These images of SST Level 4 are optimally interpolated on a global 0.01-degree grid using wavelets as basic functions. A total of 2647 images acquired between October 2006 and December 2013 were used in this study.

Front Detection with Glider Data
A gradient-based algorithm [30] was used to identify frontal zones from conservative temperature (θ fronts) and potential density (σ θ fronts) along the glider trajectory. Along-track gradients were computed using central differences and forward/backward differences [31] at locations P 1 and P 2 . The following methodology was implemented to detect frontal structures along the glider trajectory. First, glider trajectories were partitioned into three differentiated dynamical zones [15,21,24]: coastal (0-100 km offshore), transitional (100-300 km offshore), and oceanic (>300 km from the coast). Second, for each zone, a cumulative histogram of the horizontal gradients of 10-m surface-layer conservative temperature (or potential density) was computed over the study period ( Figure 2). Then, a threshold was defined for each zone as the value with a gradient magnitude greater than or equal to 90% of the cumulative frequency [10,30,32,33]; this encompasses a relatively small proportion of the distribution and includes the highest gradient values observed ( Figure 2). Finally, any θ-gradient whose value was greater than or equal to the threshold was identified as the edge of the front. When consecutive θ-gradients exceeded the threshold, the frontal interface was detected by selecting the location of the largest gradient.
where xi (i = 1, 2, ..., N) are the longitudes of the N positions along the trajectory and a0, a1, a2, a3 are the fitted coefficients for a given trajectory. Equation (2) calculates the shortest distance between the glider position and the straight line joining P1 and P2. Once the coordinates are obtained for the f (xi) = 1, 2, 3... M points, the measured SST value closest to the line (i.e., the shortest Euclidean Distance, DE) is obtained as follows: which measures the distance to the adjusted points (xa,ya) in km. Finally, data values for each trajectory were organized in a matrix containing temperature and density values, distance, and date.

MUR Observations
Global SST images from the Multi-scale Ultra-high Resolution (MUR) data product (http://mur.jpl.nasa.gov, accessed on 30 May 2019) were used for this study to compare fronts identified using in situ glider data and fronts detected from remote satellite observations [29]. These daily, cloud-free, 1-km spatial resolution images combine measurements from various satellites (microwave and infrared) and in situ observations. These images of SST Level 4 are optimally interpolated on a global 0.01-degree grid using wavelets as basic functions. A total of 2647 images acquired between October 2006 and December 2013 were used in this study.

Front Detection with Glider Data
A gradient-based algorithm [30] was used to identify frontal zones from conservative temperature (θ fronts) and potential density (σθ fronts) along the glider trajectory. Alongtrack gradients were computed using central differences and forward/backward differences [31] at locations P1 and P2. The following methodology was implemented to detect frontal structures along the glider trajectory. First, glider trajectories were partitioned into three differentiated dynamical zones [15,21,24]: coastal (0-100 km offshore), transitional (100-300 km offshore), and oceanic (>300 km from the coast). Second, for each zone, a cumulative histogram of the horizontal gradients of 10-m surface-layer conservative temperature (or potential density) was computed over the study period ( Figure 2). Then, a threshold was defined for each zone as the value with a gradient magnitude greater than or equal to 90% of the cumulative frequency [10,30,32,33]; this encompasses a relatively small proportion of the distribution and includes the highest gradient values observed ( Figure 2). Finally, any θ-gradient whose value was greater than or equal to the threshold was identified as the edge of the front. When consecutive θ-gradients exceeded the threshold, the frontal interface was detected by selecting the location of the largest gradient.  Figure 1). An upper threshold at 90% of the cumulative frequency (vertical red dashed line) was used for frontal detection. θ gradients greater than the threshold are positive detections. The same methodology was used for along-track potential density and the MUR SST series for the three study zones.

Front Detection with MUR Data
Two front detection methods were developed with MUR data. The first (MUR1D) consisted of interpolating the SST data to the coordinates of each glider trajectory during the study period. Thus, we obtained a new SST dataset fully comparable with in situ θ. The methodology described in Section 2.2.1 was applied to obtain the frontal structure with MUR1D. For the second method (MUR2D), the Sobel gradient operator was used to estimate the magnitude of the SST gradient in the images [34] since it allows to highlight this characteristic of the ocean effectively [18,35,36]. The Sobel operator convolves two 3 × 3 matrices or kernels with the SST image as follows: where SST is the MUR data, Gx and Gy are two images that contain approximations for derivatives in horizontal and vertical directions, and * indicates the convolution operator. Then, the magnitude of the SST gradients (GM) were calculated as follows: We built masks delimiting the three zones of interest: coastal (~44 km 2 ), transitional (~49.43 km 2 ), and oceanic (~49.64 km 2 ), described in Section 2.1 and represented in Figure 3a. The daily SST gradients were computed using Equation (5) (Figure 3a) for each zone. SST gradients obtained during the study period were gathered on a cumulative histogram, and a threshold value was computed following the same procedure as for the glider θ-gradients. Pixels with SST gradient intensity above the threshold correspond to a frontal zone. Images were converted into binary data, i.e., 0 (negative detection) or 1 (positive detection). Pixels that are connected were grouped and indexed with a unique label. This labeling is a fundamental processing step to differentiate frontal structures from each other [37]. When consecutive frontal points with the same label value were found, the location of the maximum SST gradient was used to define the front location Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 19 Figure 2. Histograms (bars) and cumulative frequency (red line) of the along-track gradients of conservative temperature (θ) in the coastal zone for the whole study period along CUGN Line 90 (Figure 1). An upper threshold at 90% of the cumulative frequency (vertical red dashed line) was used for frontal detection. θ gradients greater than the threshold are positive detections. The same methodology was used for along-track potential density and the MUR SST series for the three study zones.

Front Detection with MUR Data
Two front detection methods were developed with MUR data. The first (MUR1D) consisted of interpolating the SST data to the coordinates of each glider trajectory during the study period. Thus, we obtained a new SST dataset fully comparable with in situ θ.
The methodology described in Section 2.2.1 was applied to obtain the frontal structure with MUR1D. For the second method (MUR2D), the Sobel gradient operator was used to estimate the magnitude of the SST gradient in the images [34] since it allows to highlight this characteristic of the ocean effectively [18,35,36]. The Sobel operator convolves two 3 × 3 matrices or kernels with the SST image as follows: where SST is the MUR data, Gx and Gy are two images that contain approximations for derivatives in horizontal and vertical directions, and * indicates the convolution operator. Then, the magnitude of the SST gradients (GM) were calculated as follows: We built masks delimiting the three zones of interest: coastal (~44 km 2 ), transitional (~49.43 km 2 ), and oceanic (~49.64 km 2 ), described in Section 2.1 and represented in Figure  3a. The daily SST gradients were computed using Equation (5) (Figure 3a) for each zone. SST gradients obtained during the study period were gathered on a cumulative histogram, and a threshold value was computed following the same procedure as for the glider θ -gradients. Pixels with SST gradient intensity above the threshold correspond to a frontal zone. Images were converted into binary data, i.e., 0 (negative detection) or 1 (positive detection). Pixels that are connected were grouped and indexed with a unique label. This labeling is a fundamental processing step to differentiate frontal structures from each other [37]. When consecutive frontal points with the same label value were found, the location of the maximum SST gradient was used to define the front location

Front Frequency
Once the binary matrix (zeros and ones) of the frontal detection was constructed on each grid point of the transect (3 km resolution), the Front Frequency (FF) was also computed [15] for both CUGN (θ fronts and σ θ fronts) and MUR (MUR1D and MUR2D) datasets. Thus, FF can be interpreted as the probability of detecting a front along the (70) glider trajectories, expressed as follows: where N F is the number of fronts detected, and N T is the total number of trajectories. FF values obtained for each point during the study period were averaged per zone everỹ 20 km from the coast to the open ocean. Bimonthly climatologies were also constructed with the FF obtained for the in situ and satellite data sets.

Thermohaline Compensation in Density and Temperature Fronts
Thermohaline density-compensated fronts can be detected following the methodology introduced by Mauzole et al. [38]: where C i is the compensation index, α is the thermal expansion coefficient, ∇SST is the θ gradient, ∇ρ is the σ θ gradient, and r (10 3 ) is a constant introduced to obtain values of C i = 1, corresponding to the case of non-compensation. Equation (6) was calculated in the mixed layer [39] at~3 km before and after each temperature (θ) front detected. There are three cases concerning the values obtained from C i : 1. C i = 1: no thermohaline compensation; 2.

Similarity between Potential Temperature and SST
A comparison of the SST values derived from the two datasets (θ and SST) along the 70 trajectories showed a high linear Pearson's correlation coefficient (R = 0.95) and a root mean squared error (RMSE) of 0.77 • C ( Figure 4). The agreement between in situ glider and remote satellite datasets suggests that both estimates of SST across the study area are robust and contain a similar variability. The linear regression shows an intercept at 1.34 • C and R 2 = 0.903. The constant difference of 1.34 • C between the two datasets is related to (i) the differences between observational platforms and sensors used to measure the upper ocean temperature (i.e., thermistor vs. radiometer, satellite vs. in situ, etc.) and (ii) the differences between the depth of measurements.
The mean θ and SST, computed from the glider and MUR datasets as a function of offshore distance, also show a good agreement ( Figure 5). The difference in RMSE between the two mean temperatures is 0.81 • C, with the largest differences (0.86 • C) occurring near the coast (0-50 km) and the smallest (0.2 • C) occurring offshore (near 500 km). The averaged SST is always higher than the averaged computed from in situ glider data, mainly due to differences in the depth of the measurements used to compute the mean. Both the averaged θ and SST show a similar trend along the CUGN Line 90. Near the coast (0-50 km), the mean temperature decreases smoothly at 0.006 • C km −1 . Between 50 and 150 km, the mean temperature decreases rapidly (0.016 • C km −1 ), reaching a minimum near 150-200 km offshore that matches the location of the core of the California Current, which shows its maximum seasonal variability in this zone [21]. Beyond 200 km, mean temperatures increase towards the open sea but at a smaller rate (~0.003 • C km −1 ) than the one recorded between the coastal and transition zones. Monthly means showed a high correlation (R 2 = 0.865, p < 0.05, RMSE = 0.485 • C) and similar seasonal variations ( Figure A1), with a more pronounced decrease in θ and SST in spring near 150-200 km offshore. The mean θ and SST, computed from the glider and MUR datasets as a function of offshore distance, also show a good agreement ( Figure 5). The difference in RMSE between the two mean temperatures is 0.81 °C, with the largest differences (0.86 °C) occurring near the coast (0-50 km) and the smallest (0.2 °C) occurring offshore (near 500 km). The averaged SST is always higher than the averaged computed from in situ glider data, mainly due to differences in the depth of the measurements used to compute the mean.
Both the averaged θ and SST show a similar trend along the CUGN Line 90. Near the coast (0-50 km), the mean temperature decreases smoothly at 0.006 °C km −1 . Between 50 and 150 km, the mean temperature decreases rapidly (0.016 °C km −1 ), reaching a minimum near 150-200 km offshore that matches the location of the core of the California Current, which shows its maximum seasonal variability in this zone [21]. Beyond 200 km, mean temperatures increase towards the open sea but at a smaller rate (~0.003 °C km −1 ) than the one recorded between the coastal and transition zones. Monthly means showed a high correlation (R 2 = 0.865, p < 0.05, RMSE = 0.485 °C) and similar seasonal variations ( Figure   A1), with a more pronounced decrease in θ and SST in spring near 150-200 km offshore.   . Scatter plot between conservative temperature values obtained from glider at the shallowest bin (z = −10 m) and MUR SST for the whole study period. The red line is the best-fit straight line (GLIDER SST = 1.341 + 0.891 * MUR SST) for the data. N is the total number of data included in this analysis, and R 2 is the squared Pearson's correlation coefficient.
The mean θ and SST, computed from the glider and MUR datasets as a function of offshore distance, also show a good agreement ( Figure 5). The difference in RMSE between the two mean temperatures is 0.81 °C, with the largest differences (0.86 °C) occurring near the coast (0-50 km) and the smallest (0.2 °C) occurring offshore (near 500 km). The averaged SST is always higher than the averaged computed from in situ glider data, mainly due to differences in the depth of the measurements used to compute the mean.
Both the averaged θ and SST show a similar trend along the CUGN Line 90. Near the coast (0-50 km), the mean temperature decreases smoothly at 0.006 °C km −1 . Between 50 and 150 km, the mean temperature decreases rapidly (0.016 °C km −1 ), reaching a minimum near 150-200 km offshore that matches the location of the core of the California Current, which shows its maximum seasonal variability in this zone [21]. Beyond 200 km, mean temperatures increase towards the open sea but at a smaller rate (~0.003 °C km −1 ) than the one recorded between the coastal and transition zones. Monthly means showed a high correlation (R 2 = 0.865, p < 0.05, RMSE = 0.485 °C) and similar seasonal variations ( Figure   A1), with a more pronounced decrease in θ and SST in spring near 150-200 km offshore.

Temperature Gradients
Our front detection method is based on temperature (θ and SST) and potential density gradients (see Section 2). Capturing most of the gradient variability is essential for accurate front detection. Figure 6 shows the comparison between average θ, SST, and σ θ gradients for each of the study areas along Line 90 using both glider and MUR data. Despite exhibiting large standard deviations, there were significant differences (p < 0.05, Kruskal-Wallis) between the averaged θ and SST gradients obtained from the two datasets for the three zones along each trajectory. In the coastal zone, the average θ gradient from the glider data is greater (~0.05 • C km −1 , Figure 6a) than those obtained with the MUR1D (~0.02 • C km −1 , Figure 6c) and MUR2D (~0.03 • C km −1 , Figure 6d) data. This is consistent with the higher threshold obtained for the glider data (see Table 1). Using the glider data, three different zones are observed, identified in both the θ and σ θ gradients, being more intense and with greater variability in the coastal zone (0 to 100 km), where the CCS upwelling system is prone to baroclinic instabilities, driving strong mesoscale and sub-mesoscale activity. Note that the offshore decrease in the mean gradient intensity is consistent with the decrease in the threshold (see Table 1). Using MUR1D and MUR2D (Figure 6c,d), the greatest variability is found toward the transition zone (100 to 300 km).
Kruskal-Wallis) between the averaged θ and SST gradients obtained from the two datasets for the three zones along each trajectory. In the coastal zone, the average θ gradient from the glider data is greater (~0.05 °C km −1 , Figure 6a) than those obtained with the MUR1D (~0.02 °C km −1 , Figure 6c) and MUR2D (~0.03 °C km −1 , Figure 6d) data. This is consistent with the higher threshold obtained for the glider data (see Table 1). Using the glider data, three different zones are observed, identified in both the θ and σθ gradients, being more intense and with greater variability in the coastal zone (0 to 100 km), where the CCS upwelling system is prone to baroclinic instabilities, driving strong mesoscale and sub-mesoscale activity. Note that the offshore decrease in the mean gradient intensity is consistent with the decrease in the threshold (see Table 1). Using MUR1D and MUR2D (Figure 6c,d), the greatest variability is found toward the transition zone (100 to 300 km).

Front Detection
The persistence of fronts, or FF, is important for biogeochemical processes, as any front that lingers for a significant period of time may be associated with upper ocean convergence and accumulation of high biological productivity [40]. Figure 7 shows the typical distribution of FF along Line 90 as determined from CUGN (θ fronts and σθ fronts) and MUR SST (MUR1D and MUR2D) datasets. With the glider data, θ fronts and σθ fronts depict FF peaks in the coastal zone, located approximately at 100 km, with FF of 8.81 and

Front Detection
The persistence of fronts, or FF, is important for biogeochemical processes, as any front that lingers for a significant period of time may be associated with upper ocean convergence and accumulation of high biological productivity [40]. Figure 7 shows the typical distribution of FF along Line 90 as determined from CUGN (θ fronts and σ θ fronts) and MUR SST (MUR1D and MUR2D) datasets. With the glider data, θ fronts and σ θ fronts depict FF peaks in the coastal zone, located approximately at 100 km, with FF of 8.81 and 8.33, respectively. SST peaks for MUR1D and MUR2D coincide in the first 20 km near the coast, with high variability in the case of MUR1D (4.58 standard deviations). In the transition zone, the peaks of FF in θ fronts, σ θ fronts, and MUR2D are approximately 134 km from the coast. MUR2D matches the glider data FF peaks but with high variability (3.5 standard deviations). In the case of MUR1D, the peak of FF is located at 153 km with a standard deviation of 3.18. There were two FF peaks in the θ front and σ θ front data in the oceanic zone, the first at approximately 305 km, which coincides with MUR1D, and the second near 400 km. For MUR2D, a peak of FF was found at 381 km in the oceanic zone, with a mean and standard deviation of 5 and 3.58, respectively. sition zone, the peaks of FF in θ fronts, σθ fronts, and MUR2D are approximately 134 km from the coast. MUR2D matches the glider data FF peaks but with high variability (3.5 standard deviations). In the case of MUR1D, the peak of FF is located at 153 km with a standard deviation of 3.18. There were two FF peaks in the θ front and σθ front data in the oceanic zone, the first at approximately 305 km, which coincides with MUR1D, and the second near 400 km. For MUR2D, a peak of FF was found at 381 km in the oceanic zone, with a mean and standard deviation of 5 and 3.58, respectively. To compare the frontal detections obtained from the two datasets, a bimonthly climatology of FF for the study period was constructed to compare the frontal detections obtained from the two datasets ( Figure 8). Differences in FF are noticeable, with FF peaks between July and August in the estimates based on the CUGN dataset (θ fronts and σθ fronts), and between November and December in those based on the MUR dataset (MUR1D SST and MUR2D SST). The CUGN climatology of FF suggests that the front frequency follows a seasonal cycle with maxima occurring in summer. In contrast, using MUR climatology of FF shows that FF increases towards autumn. To compare the frontal detections obtained from the two datasets, a bimonthly climatology of FF for the study period was constructed to compare the frontal detections obtained from the two datasets ( Figure 8). Differences in FF are noticeable, with FF peaks between July and August in the estimates based on the CUGN dataset (θ fronts and σ θ fronts), and between November and December in those based on the MUR dataset (MUR1D SST and MUR2D SST). The CUGN climatology of FF suggests that the front frequency follows a seasonal cycle with maxima occurring in summer. In contrast, using MUR climatology of FF shows that FF increases towards autumn. In total, 703 θ fronts and 653 σθ fronts were obtained from the glider data, with ~62% of θ fronts matching σθ fronts. The oceanic region, characterized by a deeper mixed layer (~42 m) compared to the coastal (~15 m) and transition (~26 m) zones, has, on average, a higher compensation index and a high percentage of fronts that are partially compensated ( Table 2). On the other hand, the coastal zone presents the highest number of non-thermohaline-compensated θ fronts. In general, temperature fronts that coincide with density fronts do not exhibit thermohaline compensation. In total, 703 θ fronts and 653 σ θ fronts were obtained from the glider data, with~62% of θ fronts matching σ θ fronts. The oceanic region, characterized by a deeper mixed layer (~42 m) compared to the coastal (~15 m) and transition (~26 m) zones, has, on average, a higher compensation index and a high percentage of fronts that are partially compensated (Table 2). On the other hand, the coastal zone presents the highest number of non-thermohaline-compensated θ fronts. In general, temperature fronts that coincide with density fronts do not exhibit thermohaline compensation. Table 2. Averaged thermohaline density compensation index and mixed layer depth over all fronts detected using the CUGN datasets of potential density and conservative temperature for the three study zones.

Zone
Average

Interpretation of Front Detection in the CCS
We proposed and tested, for the first time, a front detection method that applies SST image processing techniques to more than seven years of conservative temperature and potential density data collected using underwater gliders off the Southern California coast. Our method analyzes cumulative histograms of temperature gradients to define a threshold used to discriminate between frontal and non-frontal observations. The threshold is computed for three different zones (coastal, transitional, and oceanic), and its values reveal marked differences among the zones, both when using CUGN (temperature and density) and MUR2D (SST) datasets. This coincides with the work of Lynn and Simpson [21], who found different physical characteristics between coastal, transitional, and oceanic zones. The largest gradients (and thresholds) are obtained in the coastal zone from the CUGN dataset. This finding suggests that frontogenesis is enhanced near the coast, leading to the intensification of horizontal gradients and driving stronger fronts.
Some peaks in the FF distribution computed from the CUGN and MUR datasets were in agreement in some areas. These areas are highly influenced by strong physical and biological interactions [21]. The peaks in FF in the coastal zone evidenced using the CUGN dataset but barely seen using the MUR dataset likely correspond to periods of intense upwelling resulting from persistent strong equatorward wind stress and a negative wind stress curl [18,41]. The FF peaks and troughs in the transitional zone are largely driven by the interaction between mesoscale energetic meanders and eddies that interact with the core of the California Current [21]. Given the large spatial and temporal variability in the transitional zone imprinted by the mesoscale, the FF peaks observed in both CUGN and MUR datasets suggest the persistence of the Ensenada Front [15]. This front is characterized by smooth surface horizontal gradients in physical, chemical, and biological oceanographic variables and is formed by a complex set of currents, some of which feed eddies [42].
The seasonal maxima in FF obtained from the CUGN datasets are consistent with stronger (weaker) gradients of dynamic ocean height occurring in spring and summer (winter) [21]. Similar consistencies were found between the bimonthly climatology of FF obtained using the CUGN datasets and the results obtained by Castelao et al. [18], who detected ocean fronts in the CCS analyzing four years of SST imagery and found the highest FF in summer (July-August) and the lowest in winter (January-February).

Differences in Front Detection with CUGN and MUR Datasets
Although the results obtained with CUGN and MUR are consistent to some extent (see Section 4.1), large discrepancies are also evidenced. In particular, the MUR datasets provide weaker temperature gradients relative to the CUGN datasets, which leads to a weaker coastal FF. This finding contrasts with Oram et al. [10], who found frontal SST gradients between 0.15 and 0.3 • C km −1 in this region. In our study, such values were attained only when in situ temperature glider data were used. In addition, the bimonthly FF obtained with MUR datasets (MUR1D and MUR2D) showed the largest FF towards the autumn, which is not consistent with the work of Castelao et al. [18], although both analyses are based on satellite SST images. This discrepancy supports the idea that the actual front dynamics is captured by CUGN but not by MUR datasets. To explain the discrepancies obtained with the two datasets, the differences in these methodological approaches are discussed below.

Temperature Sampling Depth
The MUR temperatures have been shown to present a warm bias compared to the glider data. This is because satellite SST data provide temperatures within the first meter of the water column (skin layer), which are usually warmer than the surface layer. In addition, horizontal SST gradients within the skin layer are highly sensitive to surface thermohaline-modifying processes such as surface heat fluxes (cooling/heating) and wind and wave mixing that mask front detection. On the contrary, the glider data overcome this limitation by averaging the temperature over the uppermost 10 m of the water column.

Gradient Computation
The FF estimated using MUR2D is computed from 2D horizontal gradients over the whole North Pacific Ocean off southern California. This differs from glider data, which only provide along-track θ and σ θ gradients, thus severely limiting this method. Being approximately perpendicular to the coast, the glider data along the CUGN Line 90 captures only cross-shore gradients, while the alongshore gradients are dramatically underestimated. The MUR data were also processed to conserve only the along-track SST gradients (MUR1D) to eliminate this difference. Table 3 shows a statistical comparison (Kruskal-Wallis test) of 20-km mean FF values computed using the CUGN and MUR datasets for each of the three zones. No differences were evident in the transitional and coastal zones, and the null hypothesis could not be rejected. However, there was a significant difference (p < 0.05) between MUR2D and the FF values obtained using the CUGN datasets in the oceanic zone. In this zone, there are significant differences (p < 0.05) between MUR1D and MUR2D, while the results using the MUR1D and CUGN datasets are similar in the oceanic zone. These differences between MUR and CUGN are related to the unsolved alongshore SST gradients in the glider data. Notice, however, that the FF estimated using the MUR1D datasets is not as large as the FF computed from the CUGN datasets. This is related to the weaker temperature gradient in the MUR datasets. Table 3. Comparison (Kruskal-Wallis test) of 20-km mean front frequency (FF) values estimated using the CUGN and MUR datasets for each of the three study zones. X 2 indicates the value of the Pearson's Chi-square test, while p-Value represents the level of statistical significance. Asterisks denote significant differences (p < 0.05). θ and σ θ fronts refer to fronts identified using conservative temperature ( • C) and potential density anomaly (kg m −3 ). MUR1D and MUR2D indicate the two front detection methods used here, the interpolation of SST data to the geographical coordinates of each glider trajectory to detect fronts, and the use of the Sobel gradient operator to estimate SST gradients, respectively.

Cloud Cover in MUR Data
The weaker temperature gradients in MUR are consistent with the smoothing (interpolation) procedure used to obtain 2D SST images onto a regular grid and to reconstruct the SST images during periods of severe cloud cover in MUR [43]. Indeed, the MUR2D product is a blend of satellite and in situ data interpolated onto a 1 km × 1 km grid through optimal interpolation in space or time to fill gaps where no data are available. In the presence of clouds, only data from microwave radiometers (which have low spatial resolution) and in situ data (if available for that day) are optimally interpolated, leading to a smoothed SST field compared to areas with no cloud cover. This can be observed when comparing the SST gradients reported in McClatchie et al. [23], detected with the Sobel gradient operator in level-2 images with MUR2D SST gradients ( Figure A3). On the contrary, gliders provide in situ, along-track high-resolution data regardless of cloud conditions. The cloud cover limitation in MUR becomes critical when selecting detection algorithms since, in addition to detecting extreme gradients, the algorithm should not be overly sensitive to the noise caused by missing data to avoid false detections.
To explain the difference in the FF of climatological cycles from the MUR and glider data, a 20-year climatology of the total cloud cover distribution in the North Pacific Ocean off southern California ( Figure A2) suggests that the offshore (oceanic) zone of CUGN Line 90 is largely covered by clouds most of the year, especially in July to August. In contrast, the coastal zone has the largest cloud cover during January to February. This might explain the low FF computed using the MUR datasets in summer (winter) in the oceanic (coastal) zones (Figures 7 and 9).

Advantages and Disadvantages of Detecting Fronts from MUR and Glider Datasets
Satellite SST data have worldwide coverage and thus offer the advantage of being available with a large spatio-temporal range, avoiding the need to conduct in situ expensive (time and money) campaigns. In addition, they allow computing 2D SST gradients. However, their resolution is considerably reduced in cloudy areas where infrared radiation cannot be measured. In these areas, the spatial smoothing generated by the MUR optimal interpolation reduces the value of spatial temperature gradients. Glider observations, by contrast, provide a high-resolution but challenging dataset that convolves spatial, temporal, or spatio-temporal variability as the glider moves along the CalCOFI Line 90. In addition, glider data provide vertical profiles. Thus, they allow cross-track front detection in deeper layers and the study of the vertical structure of fronts ( Figure 10). Our method presented in this paper can be easily applied to other depths.

Advantages and Disadvantages of Detecting Fronts from MUR and Glider Datasets
Satellite SST data have worldwide coverage and thus offer the advantage of being available with a large spatio-temporal range, avoiding the need to conduct in situ expensive (time and money) campaigns. In addition, they allow computing 2D SST gradients. However, their resolution is considerably reduced in cloudy areas where infrared radiation cannot be measured. In these areas, the spatial smoothing generated by the MUR optimal interpolation reduces the value of spatial temperature gradients. Glider observations, by contrast, provide a high-resolution but challenging dataset that convolves spatial, temporal, or spatio-temporal variability as the glider moves along the CalCOFI Line 90. In addition, glider data provide vertical profiles. Thus, they allow cross-track front detection in deeper layers and the study of the vertical structure of fronts ( Figure 10). Our method presented in this paper can be easily applied to other depths. tion cannot be measured. In these areas, the spatial smoothing generated by the MUR optimal interpolation reduces the value of spatial temperature gradients. Glider observations, by contrast, provide a high-resolution but challenging dataset that convolves spatial, temporal, or spatio-temporal variability as the glider moves along the CalCOFI Line 90. In addition, glider data provide vertical profiles. Thus, they allow cross-track front detection in deeper layers and the study of the vertical structure of fronts ( Figure 10). Our method presented in this paper can be easily applied to other depths.

Sensitivity of the Results to the Front Definition
Our results rely on the definition of a frontal zone. In our method, a front is defined as gradients stronger than the 90th percentile in the study area. The value of this threshold thus depends on the area studied and the time covered by the dataset used. Longer records, such as the one used in this study, provide a robust detection of frontal structures since they capture a large portion of the temporal variability in temperature (and density) gradients at a wide range of scales. It should be noted that the threshold is objectively selected as the value at the inflection point on the cumulative function of temperature gradients [44]. This was evidenced by Oram et al. [10], who used the Canny algorithm [8] to detect fronts in SST images and recommended defining upper and lower thresholds corresponding to 90% (as in our study) and 80%, respectively, in the cumulative function of SST gradients. Computing the threshold for the area studied is necessary because the SST gradient intensity is highly variable worldwide. Furthermore, the use of a fixed threshold value would imply an under-or overestimation of the number of frontal features detected if dynamical conditions change.
The use of the thresholds defined in this study plays an important role in the detection of fronts because it can define regions that are affected by various oceanographic processes. Saraceno et al. [30] applied different thresholds for six regions of the Southwestern Atlantic Ocean using SST data, where frontogenesis associated with larger-scale deformation increases the value of the SST threshold and the variability of the FF.

Differences between Temperature and Density Fronts: Thermohaline Compensation
The upper-ocean mixed layer is often characterized by thermohaline compensation. The contributions of temperature and salinity to density tend to (partially) cancel each other out [45]. Thermohaline compensation plays a central role in detecting ocean fronts from density and may explain the differences between fronts detected from SST and density data. Fronts detected with temperature gradients agree with fronts detected with density gradients when there is no thermohaline compensation [38].
In this work, it was indeed observed that most of the thermal fronts detected have no thermohaline compensation index coinciding with density fronts (Figure 11a,b). When surface temperature gradients are weak but with no thermohaline compensation, the thermal front may go undetected (uppermost value below the threshold), although the density front may sometimes be detected (Figure 11c,d). A deepening of the mixing layer has been shown to strengthen the thermohaline compensation [46]. The mixed layer is deeper in the oceanic zone in our study region, where the weakest temperature gradients are found relative to the transition and coastal zones. The coastal zone shows a shallow mixing layer where temperature contributes the most to density gradients [46], so that most of the thermal fronts detected in this region are not thermohaline compensated and are also observed in density.

Conclusions
We used an automated method to detect ocean fronts from continuous glider measurements and satellite imagery. The method is based on determining a threshold value for the horizontal gradients of temperature and density from a cumulative climatological function. This threshold has been computed for the three zones studied: coastal, transition, and oceanic. Three different threshold values have been obtained, depending on the particular dynamics of each zone.
FF values computed using CUGN glider data and high-resolution MUR SST images were compared, and the spatial and temporal variability along the CUGN Line 90 was analyzed. There was moderate agreement in the identification of frontal peaks in the FF distribution concerning the distance from the coast in the three regions. However, strong differences were found in the bimonthly FF probability density functions of the computed gradients, with high (low) FF in spring-summer for the glider (MUR) data. Such differences are associated with the smoothing of the MUR data resulting from the objective interpolation, the difference in the depth of measurements, and the orientation of fronts related to the glider-track direction.
The most relevant conclusions regarding front frequency in the three zones studied are as follows: • Front Frequency (and horizontal gradients) increases significantly towards the coastal zone and in summer; • The bimonthly climatology of FF using CUGN glider data shows a seasonal cycle, i.e., high in July-August and low in January-February, while the MUR dataset shows high FF values in November-December and low in May-June; • The oceanic zone, distant from the CC core, exhibits the weakest horizontal gradients of the three zones, and thus fewer frontal structures are detected in this area;

Conclusions
We used an automated method to detect ocean fronts from continuous glider measurements and satellite imagery. The method is based on determining a threshold value for the horizontal gradients of temperature and density from a cumulative climatological function. This threshold has been computed for the three zones studied: coastal, transition, and oceanic. Three different threshold values have been obtained, depending on the particular dynamics of each zone.
FF values computed using CUGN glider data and high-resolution MUR SST images were compared, and the spatial and temporal variability along the CUGN Line 90 was analyzed. There was moderate agreement in the identification of frontal peaks in the FF distribution concerning the distance from the coast in the three regions. However, strong differences were found in the bimonthly FF probability density functions of the computed gradients, with high (low) FF in spring-summer for the glider (MUR) data. Such differences are associated with the smoothing of the MUR data resulting from the objective interpolation, the difference in the depth of measurements, and the orientation of fronts related to the glider-track direction.
The most relevant conclusions regarding front frequency in the three zones studied are as follows: • Front Frequency (and horizontal gradients) increases significantly towards the coastal zone and in summer; • The bimonthly climatology of FF using CUGN glider data shows a seasonal cycle, i.e., high in July-August and low in January-February, while the MUR dataset shows high FF values in November-December and low in May-June; • The oceanic zone, distant from the CC core, exhibits the weakest horizontal gradients of the three zones, and thus fewer frontal structures are detected in this area; • Thermohaline-compensated fronts are more abundant towards the oceanic zone, although most fronts are detected using both temperature or density criteria, indicating that the contribution of temperature to density in this region is the most important thermohaline contribution. Figure A2. Bimonthly climatology of total cloud cover from ERA5 for the study area between 1 January 1994 and 12 December 2013.